An Approximation Algorithm for Computing 
Shortest Paths in Weighted 3-d Domains* 

Lyudmil Aleksandrov^ Hristo Djidjev * Anil Maheshwari 

Jorg-Riidiger Sack ^ 

O 

<N : January 12, 2013 

in 



u 



Abstract 

We present the first polynomial time approximation algorithm for computing short- 
est paths in weighted three-dimensional domains. Given a polyhedral domain D, con- 
sisting of n tetrahedra with positive weights, and a real number e 6 (0, 1), our algorithm 
constructs paths in T> from a fixed source vertex to all vertices of T>, whose costs are 
at most 1 + e times the costs of (weighted) shortest paths, in 0(C(T>)-j^ log j log 3 ^) 
time, where C(T>) is a geometric parameter related to the aspect ratios of tetrahedra. 
J> , The efficiency of the proposed algorithm is based on an in-depth study of the 

local behavior of geodesic paths and additive Voronoi diagrams in weighted three- 
dimensional domains, which are of independent interest. The paper extends the results 
C*~) . of Aleksandrov et al. [1] to three dimensions. 

1 Introduction 
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1.1 Motivation 

The computation of shortest paths is a key problem arising in a number of diverse applica- 
tion areas including geographic information systems, robotics, computer graphics, computer- 
aided design, medical computing and others. This has motivated the study and subsequent 
design of efficient algorithms for solving shortest path problems in different settings based on 
the geometric nature of the problem domain (e.g., two-dimensional (2-d), three-dimensional 
(3-d), surfaces, presence/absence of obstacles) and the cost function/metric (e.g., Euclidean, 
L p , link distance, weighted/unweighted, multi-criteria). In addition to its driver - the appli- 
cations - the field has provided, and continues to do so, exciting challenges from a theoretical 
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perspective. As a result, shortest path problems have become fundamental problems in areas 
of Computer Science such as Computational Geometry and Algorithmic Graph Theory. 

The standard 3-d Euclidean shortest path problem of computing a shortest path between 
pair of points avoiding a set of polyhedral obstacles, denoted as the ESP3D problem, is known 
to be iVP-hard even when the obstacles are parallel triangles in the space. It is not difficult 
to see that the number of combinatorially distinct shortest paths from a source point to a 
destination point may be exponential in the input size. Canny and Reif [5] used this to 
establish the iVP-hardness of the ESP3D problem, for any L p metric, p > 1. In addition to 
this combinatorial hardness result, Bajaj [7] has provided an algebraic hardness argument 
that an exponential number of bits may be required. More recently, Mitchell and Sharir [22] 
gave iVP-completeness proofs for the problem of computing Euclidean shortest paths among 
sets of stacked axis-aligned rectangles, and computing Li-shortest paths among disjoint balls. 
Given the AP-hardness of the ESP3D problem, work has focused on exploiting the geometric 
structure of the obstacles and/or on providing approximation algorithms. We will mention 
some of these approaches in Section 11.41 

In many applications, the Euclidean metric does not capture adequately the nature of 
the problem, for instance when the problem domain is not homogeneous. This motivates 
the weighted versions of the shortest path problem. For example in the 2-d case, consider 
triangulated regions where each triangle represents a particular terrain type such as water, 
rock, or forest. Here different weights capture the cost of traveling a Euclidean unit-length 
through each face. Incorporating weights makes the solution more difficult to obtain even 
in 2-d, but it does provide more realistic answers. It is known that light and other types 
of waves (e.g., seismic and sonic) travel along the shortest paths in heterogeneous media. 
Hence, algorithms solving the weighted shortest path problem (WSP3D) can be used for 
modeling wavefront propagation in such media. In the 3-d, a number of applications are 
non-homogeneous in nature and can be expressed using the weighted model. Next, we list 
some of such potential applications. 

• In geology, seismic refraction and reflection methods are used based on measurements of 
the travel time of seismic waves refracted at the interfaces between subsurface layers of 
different densities. As such waves propagate along shortest paths and weighted shortest 
path algorithms may be used to produce more accurate and more efficient estimation 
of subsurface layer characteristics, e.g., the amount of oil contained in the subsurface 
[14] . Another related application is the assessment of garbage dumps' health. When 
a new garbage dump is built, sensors are placed at the bottom, and when the garbage 
dump starts to fill, waves from the top passing through the garbage to these sensors 
are used in order to determine the decomposition rate of the garbage [14"] . 

• Computation of 3-d shortest path have also been used to compute fastest routes for 
aircrafts between designated origin and destination points while avoiding hazardous, 
time- varying weather systems. Krozel et al. [TTJ investigate synthesizing weather avoid- 
ance routes in the transition airspace. Our weighted 3-d region model can be used to 
generalize that approach: instead of totally avoiding undesirable regions, one can as- 
sign penalty weights to them and then search for routes that minimize travel through 
such regions, while also avoiding unnecessarily long detours. 
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• In medical applications simulation of sonic wavefront propagation is used when per- 
forming imaging methods as photoacoustic tomography or ultrasound imaging through 
heterogeneous tissue [121 [21]. In radiation therapy, domain features include densities 
of tissue, bone, organs, cavities, or risk to radiation exposure, and optimal radiation 
treatment planning takes this non-homogeneity into consideration. 

• The problem of time-optimum movement planning in 2-d and 3-d for a point robot that 
has bounded control velocity through a set of n polygonal regions of given translational 
flow velocities has been studied by Reif and Sun [21]. They state that this intriguing 
geometric problem has immediate applications to macro-scale motion planning for 
ships, submarines, and airplanes in the presence of significant flows of water or air. 
Also, it is a central motion planning problem for many of the meso-scale and micro-scale 
robots that have environments with significant flows that can affect their movement. 
They establish the computational hardness for the 3-d version of this problem by 
showing the problem to be PSPACE hard. They give a decision algorithm for the 
2-d flow path problem, which has very high computational complexity, and they also 
design an efficient approximation algorithm with bounded error. The determination of 
the exact computational complexity of the 3-d flow path problem is posed as an open 
problem. Although, our weighted 3-d model does not apply directly to this setting, 
it can be used to construct initial approximations by assigning appropriate weights 
depending on the velocity and direction of the flows in different regions. In addition, 
the discretization scheme and the algorithmic techniques developed here can be used 
for solving the 3-d flow path problem. 

1.2 Problem formulation 

In this paper, we consider the following problem. Let Dbea connected 3-d domain consisting 
of n tetrahedra with a positive real weight associated to each of them. The 3-d weighted 
shortest path problem (WSP3D) is to compute minimum-cost paths in D from a fixed source 
vertex to all vertices of T>. The cost of a path in T> is defined as the weighted sum of the 
Euclidean lengths of the sub-paths within each crossed tetrahedron. We will describe and 
analyze an approximation algorithm for this problem that, for any real number e G (0, 1), 
computes paths whose costs are at most 1 + e times greater than the costs of the minimum 
cost paths. In Section [21 we describe our model in detail. 

Note that the WSP3D problem can be viewed as a generalization of the ESP3D problem. 
Namely, given an instance of the ESP3D problem, one can find a large enough cube containing 
all the obstacles, tetrahedralize the free-space (i.e., exterior of the obstacles, but in the 
interior of the cube) and set equal weights to the resulting tetrahedra obtaining an instance 
of the WSP3D problem. 

1.3 Challenges 

A key difference between Euclidean shortest path computation in 2-d and 3-d weighted 
domain is the A^P-hardness already mentioned. Underlying this is the fact that, unlike in 2- 
d, Euclidean 3-d shortest paths are not discrete. Specifically, in 2-d, the edges of a shortest 
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path (e.g., Euclidean shortest paths among obstacles in the plane) are edges of a graph, 
namely, the visibility graph of the obstacles including the source and the destination points. 
In contrast, in polyhedral 3-d domains, the bending points of shortest paths on obstacles 
may lie in the interior of the obstacles' edges. Moreover, in weighted 3-d settings, bending 
points may even belong to the interior of the faces. 

Furthermore, even in the case of weighted triangulated planar domains, the (weighted) 
shortest path may cross each of the n cells Q(n) times and may be composed of 0(n 2 ) seg- 
ments in total. Not only is the path complexity higher, but the computation of weighted 
shortest paths in 2-d turns out to be substantially more involved than in the Euclidean set- 
ting. In fact, there is not even an exact algorithm known, and the first approximation 
algorithm due to [21] had an 0(n 8 log(-)) time bound, where n is the number of triangles 
in the subdivision. This problem has been actively researched since then, and currently the 
best known algorithm for the weighted region problem on planar subdivisions (as well as on 
polyhedral surfaces) runs in 0(^ log ^ log ^) time [3]. (Also, see [I] for a detailed literature 
review for the planar case.) 

One of the classical tools of Computational Geometry is the Voronoi Diagram. This struc- 
ture finds numerous applications (see e.g., [6]). It is also a key ingredient in several efficient 
shortest path algorithms. Researchers have studied these diagrams under several metrics 
(including Euclidean, Manhattan, weighted, additive, convex, abstract) and for different 
types of objects (including points, lines, curves, polygons), but somehow the computation of 
these diagrams in media with different densities (i.e., the refractive media) remained elusive. 
One of the main ingredients in solving the problem studied here is to compute (partial) 
additive Voronoi diagrams of points in refractive media. The generic techniques of Klein 
[T5| [16] and Le [19] do not apply in this the bisecting surfaces do not satisfy the 

required conditions. In this paper, we make an important step towards the understanding 
and computation of these diagrams. 



1.4 Previous related work 

By now, shortest path problems in 2-d are fairly well understood. Efficient algorithms have 
been developed for many problem instances and surveys are readily available describing the 
state of the art in the field. 

In 3-d, virtually all the work has been devoted to the ESP3D problem. Papadimitriou 
[23] suggested the first polynomial time approximation scheme for that problem. It runs 
in 0(^2 (L + \og(n/e)) time, where L is the number of bits of precision in the model of 
computation. Clarkson [TTJ provided an algorithm running in 0(n 2 X(n) \og(n / e) / (e 4 ) + 
n 2 lognplog(nlogp)) time, where p is the ratio of the longest obstacle edge to the distance 

between the source and the target vertex, A(n) = a(n)°^ n ^ , and a(n) is the inverse 
Ackermann's function. 

Papadimitriou's algorithm was revised and its analysis was refined by Choi et al. [S] un- 
der the bit complexity framework. Their algorithm runs roug hly in 0( !L 4-fj,(X)) time, 
where fi(X) represents the time (or bit) complexity of multiplying X-bit integers and 
X = 0(log(-) + L). In [TU], the same authors further developed their ideas and pro- 
posed a precision-sensitive algorithm for the ESP3D problem. In [5], Asano et al. proposed 
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and studied a technique for computing approximate solutions to optimization problems and 
obtained another precision-sensitive approximation algorithm for the ESP3D problem with 
improved running time in terms of L. 

Har-Peled [13] proposed an algorithm that invokes Clarkson's algorithm as a subroutine 
0(^2 log -) times to build a data structure for answering approximate shortest path queries 
from a fixed source in O(log-) time. The data structure is constructed in roughly 0(tt) 
time. Agarwal et al. [TJ considered the ESP3D problem for the case of convex obstacles and 
proposed an approximation algorithm running in 0(n + ^ log 2 - log log k) time, where k is 
the number of obstacles. In contrast to all other algorithms discussed here, the complexity 
of this algorithm does not depend on the geometric features of the obstacles. In the same 
paper, the authors describe a data structure for answering approximate shortest path queries 
from a fixed source in logarithmic time. 

In the weighted (non-Euclidean) 3-d case no previous algorithms have been reported 
by other authors. In [3], we have announced and sketched a polynomial time approxima- 
tion scheme for WSP3D problem that runs in 0(^5 log + logn)) time. The run-time 
improves to log Mogn) when all weights are equal. This algorithm can be used to 
efficiently solve the ESP3D problem. In this paper, we apply that approach, but develop 
the required details, apply new techniques, improve the complexity bounds, and provide a 
rigorous mathematical analysis. 

1.5 Contributions of this paper 

In this paper, we make several contributions to the fields of shortest path computations and 
the analysis of weighted 3-d regions model, as listed below. 

• We provide an approximation algorithm for solving the WSP3D problem in a poly- 
hedral domain T> consisting of n weighted tetrahedra. The algorithm computes ap- 
proximate weighted shortest paths from a source vertex to all other vertices of D 
in 0(C(V)-^ log - log 3 -) time, where e E (0,1) is the user-specified approximation 
parameter and C(T>) is a geometric parameter related to the aspect ratios of tetrahe- 

The cost of the computed paths are within a factor of 1 + e of the cost of the 
corresponding shortest paths. 

As we stated above, the ESP3D problem, i.e., the unweighted version of this problem, 
is already iVP-hard even when the obstacles are parallel triangles in the space [Hj . The 
time complexity of our algorithm, which is designed for the more general weighted 
setting, compares favorably even when applied in the Euclidean setting to the existing 
approximation algorithms. 

• Our detailed analysis, especially the results on additive Voronoi diagrams derived in 
Section [2j provides valuable insights into the understanding of Voronoi diagrams in 
heterogeneous media. This may open new avenues, for example, for designing an 
algorithm to compute discretized Voronoi diagrams in such settings. 

1 See Lemma I37L1 for details on the value of C(D). 
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• Our approximation algorithms in 2-d have proven to be easily implementable and 
of practical value [IB]. Our algorithm for WSP3D presented here, in spite of being 
hard to analyze, essentially uses similar primitives, and thus has the potential to be 
implementable, practical, and applicable in different areas. 

• Our work provides further evidence that discretization is a powerful tool when solving 
shortest path-related problems in both Euclidean and weighted settings. We conjecture 
that the discretization methodology used here generalizes to any fixed dimension. 

Furthermore, our discretization scheme is independent of the source vertex and can be 
used with no changes to approximate paths from other source vertices. This feature 
makes it applicable for solving all pairs shortest paths problem and for designing data 
structures for answering shortest path queries in weighted 3-d domains. 

• The complexity of our algorithm does not depend on the weights assigned to the 
tetrahedra composing V, but it depends on their geometry. We analyze and evaluate 
that dependence in detail. Geometric dependencies arise also in Euclidean settings and 
in most of the previous papers. For example, in Clarkson [TTJ, the running time of the 
algorithm depends on the ratio of the longest edge to the distance between the source 
and the target vertex. Applying known techniques (see e.g., pQ), such dependency can 
often be removed. Here, this would be possible provided that an upper bound on the 
number of segments on weighted shortest paths in 3-d is known. However, the increase 
in the dependency on n in the time complexity that these techniques suffer from, which 
is of order Q(n 2 ), appears not to justify such an approach here. In our approach, the 
dependency on the geometry is proportional to the average of the reciprocal squared 
sinuses of the dihedral angles of the tetrahedra composing T>. Thus, when n is large, 
many tetrahedra would have to be fairly degenerate so that this average to play a major 
role. We therefore conjecture that in typical applications, the approach presented here 
would work well. 

1.6 Organization of the paper 

In Section [21 we describe the model used throughout this paper, formulate the problem, 
present some properties of shortest paths in 3-d, and derive a key result on additive Voronoi 
diagrams in refractive media. In Section [3], we describe our discretization scheme which is 
a generalization of the 2-d scheme introduced in [3]. In Section HI we construct a weighted 
graph, estimate the number of its nodes and edges and prove that shortest paths in T> can be 
approximated by paths in the graph so constructed. In Section we present our algorithm 
for the WSP3D problem. In Section [61 we conclude this paper. 

2 Problem formulation and preliminaries 
2.1 Model 

Let D be a connected polyhedral domain in 3-d Euclidean space. We assume that T> is 
partitioned into n tetrahedra 7\, . . . ,T n , such that T> = U" =1 Tj and the intersection of any 
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pair of different tetrahedra is either empty or a common element (a face, an edge, or a vertex) 
on their boundaries. We call these tetrahedra cells. A positive weight Wi is associated with 
each cell Tj representing the cost of traveling in it. The cost of traveling along a boundary 
element of a cell is the minimum of the weights of the cells incident to that boundary element. 
We consider paths (connected rectifiable curves) that stay in T>. The cost of a path ti in T> 
is defined by ||7r|| = ^ILi^iW; where denotes the Euclidean length of the intersection 
7Tj = 7r fl Tj. Boundary edges and faces are assumed to be part of the cell from which they 
inherit their weight. 

Given two distinct points u and v in D, the shortest path problem in T> is to find a 
minimum cost path ir(u, v) between u and v that stays in T>. We refer to the minimum cost 
paths as shortest paths. For a given approximation parameter e > 0, a path tt £ = ir e (u,v) is 
an ^-approximation of the shortest path it = tt(u,v), if ||7r e || < (1 +e)||7r||. Without loss of 
generality, we may assume that the points u and v are vertices of T>, since, if they are not, we 
can make them such by partitioning the cells where they belong. In this paper, we present 
an algorithm that, for a given source vertex u and an approximation parameter e G (0, 1), 
computes e-approximate shortest paths from u to all vertices of T>. 

In this setting, it is well known |2lj§ that shortest paths are simple (non self-intersecting) 
and consist of a sequence of segments, whose endpoints are on the cell boundaries. The 
intersection of a shortest path with the interior of a cell, a face, or an edge is a set of disjoint 
segments. More precisely, each segment on a shortest path is of one of the following three 
types: 

(i) cell-crossing - a segment that crosses a cell joining two points on its boundary; 

(ii) face-using - a segment lying along a face of a cell; 

(iii) edge-using - a segment along an edge of a cell. 

We define linear paths to be paths consisting of cell-crossing, face-using, and edge-using 
segments exclusively. A linear path ir(u, v) can be represented as the sequence of its segments 
{si, . . . , S[ + i} or, equivalently, as the sequence of points {a , . . . , a>i+i}, lying on the cell 
boundaries that are endpoints of these segments, i.e., Sj = (aj_i,aj), u = ao, and v = az+i. 
The points Oj that are not vertices of cells are called bending points of the path. 

The local behavior of a shortest path around a bending point a, lying in the interior of a 
face /, is fully described by the directions of the two segments of the shortest path, s~ and 
s + , that are incident to a. The direction of each of these two segments is described by a pair 
of angles, which we denote by (ip~ , 6~) and (ip + , 9 + ), respectively. The in-angle <p~ is defined 
to be the acute angle between the direction normal to / and the segment s~ . Similarly, the 
out-angle (p + is the acute angle between the normal and the segment s + . The angles 9~ and 
9 + are the acute angles between the orthogonal projections of s~ and s + with a reference 
direction in the plane containing the face /, respectively (see Figured]). 

It is well known that when it is a shortest path it is a linear path such that the angles 
(<p~ , 0~) and (tp + , 6 + ) are related by Snell's law as follows: 



2 The 2-d case was treated there, but the arguments readily apply to the 3-d model considered in this 
paper. 
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Figure 1: An illustration of the Snell's law of refraction. 



Snell's Law of Refraction: Let a be a bending point on a geodesic path it lying in the 
interior of a face f of D. Let the segments preceding and succeeding a in tx be s~ and s + , 
respectively. Let the weights of the cells containing s~ and s + be w~ and w + , respectively. 
Then s + belongs to the plane containing s~ and perpendicular to f and the angles p~ and 
p + are related by w~ sin<^~ = w + sm<p + . 

We refer to linear paths that are locally optimal (i.e., satisfy the Snell's law) as geodesic 
paths. Hence, the shortest path between pair of vertices u and v is the geodesic path of 
smallest cost joining them. In the following we discuss some of the implications of Snell's 
law on the local behavior of geodesic paths. 

Hereafter, we denote by k the ratio w + /w~ . Without loss of generality, we assume that 
w~ > w + , i.e., k < 1. Let tp* be the acute or right angle for which sinp* = k. We refer to 
this angle as the critical angle for the face /. From Snell's law, it follows that p~ < p*. The 
case where p~ = ip* deserves a special attention. In this case, p + must be a right angle and 
therefore the segment s + is a face-using segment. Furthermore, if the second endpoint a\ of 
s + is in the interior of /, then the segment following s + is inside the tetrahedron containing 
s~, and the out-angle at a\ is equal to p* (see Figured (b)). In summary, if s is a face-using 
segment, then it is preceded and followed by segments lying in the cell with bigger weight 
and their corresponding in- angle and out-angles are equal to the critical angle p*. 

In the next subsection, we study the properties of simple geodesic paths joining points 
in neighboring cells or in the same cell through a face- using segment. We define a function 
related to the cost of these geodesic paths and prove a number of properties that it possesses. 
These properties are essential to the design and the analysis of our algorithm. 

2.2 Weighted distance function 

Let F be a plane in the three-dimensional Euclidean space. We denote the two half-spaces 
defined by F by F~ and F + and assume that positive weights w~ and w + have been asso- 
ciated with them, respectively. We extend our model by assigning a weight w to F, so that 
w = min(w~,w + ) if w~ ^ w + , and < w = w + (= w~) if w~ = w + . The latter case models 
the situation where the geodesic path joins two points in the same cell through a face-using 
segment on the boundary of that cell. 
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We refer to the half-spaces F~ and F + as the lower and the upper half-space, respectively. 
Let v be a point in the lower half-space F~ at distance z~ from F and £ be a line parallel 
to .F in the upper half-space _F + at distance z + from F. Let O xyz be a Cartesian coordinate 
system such that the plane O xy coincides with F, v has coordinates (0,y, —z~), and the line 
£ is described by {£ : y = 0, z = z + }. 

We consider a point x = (x, 0, z + ) on I and denote by 7f (v , x) the geodesic path between 
v and x. In this setting, the geodesic path is unique and thus coincides with the shortest 
path. We denote the cost of this path by c(v,x), where x is the x-coordinate of x. So, for 
fixed I, c(v, x) can be viewed as a function defined for any real x. We call c the weighted 
distance function from v to I (Figure [2]). 

The local structure of the geodesic path 7f(v,x) is governed by Snell's law. In the case 
where w~ ^ w + , the shortest path between v and x consists of two segments (v,a) and 
(a, x), where the bending point a is uniquely determined by Snell's law (Figure [3] (a)). In 
the case where w~ = w + , the structure of the path ir(v, x) is as follows. It is a single segment 
(v,x), provided that the angle if between (v,x) and the direction normal to the plane F is 
smaller than or equal to the critical angle tp* defined by simp* = w/w~. Or, if <p > tp*, it 
is in the plane perpendicular to F containing v and x and consists of three segments (v, a), 
(a, ai) and (ai,x), where the acute angles between the segments (v,a) and (ai,x), and the 
direction normal to the plane F are equal to the critical angle tp*, and the segment (a, ai) is 
in F (Figure 0(b)). 

From these observations, it follows that, in all cases, weighted distance function can 
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equivalently be defined by 

c(v,x) — ||7f(v,x)|| = min (w~|wa| + u;|aai| + w; + |aix|). (1) 

a,ai£F 

In the case where w~ ^ w + , the minimum is achieved when a = a\ = (tx, (1 — r)y, 0), where 
r is the unique solution in (0, 1) of the equation 

^ 



2\ 



yjr 2 (x 2 + y 2 ) + 0") 2 ^(1 - t) 2 (x 2 + y 

where v = (0,y,z~) and x = (x,0,z + ). The latter leads to an algebraic equation of degree 
four and it is infeasible^l to evaluate c(v,x) explicitly. 

The case where w~ = w + is easier, as in that case the geodesic path is either a straight 
line, or a three segment path as described above and illustrated in Figure [3] (b) and the 
function c(v, x) has an explicit representation, which is 

w~ \J 1 x 2 + y 2 + z 2 if ' x 2 + y 2 < z/w ^ 
w(^/x 2 + y 2 — zw) if \fx 2 ~+~y 2 > z/w, 

where z = z~ + z + and w = \J [w~ /w) 2 — 1. We refer to this case as the explicit case. In 
the next lemma we state and prove some useful properties of the function c(v,x). 

Lemma 2.1 For a fixed v, the weighted distance function c(v, x) has the following properties: 

(a) It is continuous and differentiate. 

(b) It is symmetric, i.e. c(v,x) = c(v, —x). 

(c) It is strictly increasing for x > 0. 



3 Although the roots of a quartic can be expressed as a rational function of radicals over its coefficients, 
they are too complex to be analytically manipulated and used here. 
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(d) It is convex. 

(e) It has asymptotes for x — > ±00 as follows: 

(el) if w + < w~ then the asymptotes are w + (z~ cotip* ± x), 
(e2) if w~ < w + then the asymptotes are w~(z + cotip* ± x), 
(e3) in the explicit case w + = w~ > w the asymptotes are ±wx, 
where ip* is the critical angle. 

Proof: In the explicit case (w~ = w + ), all these properties follow straightforwardly from 
the explicit representation ([3]). So, we consider the case w~ 7^ w + 

From ([1]) and a 1 = a = (tx, (l—r)y, 0) it follows that c(v , x) = w~ \/t 2 (x 2 + y 2 ) + (z~) 2 + 
w + \J (1 — t) 2 (x 2 + y 2 ) + (z + ) 2 , where r is the root of the equation ([2]). The root r can 
be viewed as a function of x, which by the implicit function theorem is continuous and 
differentiable. Hence property (a) holds. 

The property (b) follows from the observation that the value of the function c(v, x) is 
determined by the distance between the projections of the points v and x on F, which is 
\/ y 2 + x 2 where y is fixed. 

To prove (c) we consider a point x' = [x! , 0, z + ) such that x' > x > and de- 
note by t' the corresponding root of (j2J). We have c(v , x') = w~ t i2 [x' 2 + y 2 ) + (z~) 2 + 
w + \J (1 — t') 2 (x' 2 + y 2 ) + (z + ) 2 . Using the fact that the function c(v,x) is defined as the 
cost of the shortest path joining v and x we have 

c(v,x) < w~^/t' 2 (x 2 + y 2 ) + (z~) 2 + w + y/(l - t') 2 (x 2 + y 2 ) + (z + ) 2 
< w-y/r' 2 (x' 2 + y 2 ) + (z~) 2 + w + ^(l - r') 2 {x' 2 + y 2 ) + (z+) 2 = c(v, x'). 

In order to prove (d), we show that for any three equidistant points ~x.\ < x < x 2 on 
£, i.e., such that 2x = X\ + x 2 , the second finite difference A 2 (c; X\, x , x 2 ) = c(v,Xi) — 
2c(v, xo) + c(v, £2) of the function c(v, x) is positive. We denote by a± and a 2 the bending 
points of the shortest paths from v to xi and x 2 , respectively. Let a' Q be the middle point 
of the segment (01,02). Then, using the definition of c(v,xo) and the convexity of the 
Euclidean distance function we obtain 2c(v,x ) < 2(w~\va' \ + u> + |a'x |) < w~(\vai\ + 
|ua 2 |) + w + (|aiX!| + I a-2 x 2 1 ) = c(v, X\) + c{y , x 2 ), which implies A 2 (c; x\, x , x 2 ) > and (d). 

Finally, we prove (e). Let us assume that w + < w~ . In this case, using 
Snell's law we observe that when x — > +00 the bending point a(x) of the shortest 
path 7f(v,x) converges to the point [z~ tany?*, y, 0) (see Figure [2]). Hence, we have 
\im x ^ +00 (w + \x — z~ tany?*) 2 + y 2 + (z + ) 2 + w~z~ /cosip* — c(f,x)) = 0. On the other 
hand 



lim (w + a/ (x — z~ tany?*) 2 + y 2 + (z + ) 2 + w z / cosip* — w + (z cot ip* + x)) 



= w + lim (\J {x — z tan</)*) 2 + y 2 + {z + ) 2 — (x — z tan ip*)) 

x—>+oo 

= w lim — = 0. 

(y/(x - z~ tany?*) 2 + y 2 + (z + ) 2 + (x — z~ tamp*) 

Combining these two limits we obtain \im x ^ +OQ (c(v , x) — w + (z~ cotip* + x)) = and thus 
(el) is valid for x — > +00. The case where x — > — 00 is symmetric. 
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In the case where w~ < w + we use Snell's law and observe that the bending point a{x) of 
the shortest path 7f (v, x) converges to x— z + tamp* , that is lim a ,_ >+00 (a(x)— x — z + tan ip*) = 0. 
Then (e2) is established analogously to (el). □ 

2.3 Refractive Additive Voronoi diagram 

Next we study Voronoi diagrams under the weighted distance metric defined above. Given 
a set S of k points v\, . . . , vu in F~ , called sites, and nonnegative real numbers C\, . . . , Cf., 
called additive weights, the additive Voronoi diagram for 5* is a subdivision of F + space into 
regions V{yi,F + ) = {x G F + | dist(x, Vi) + Ci < dist(x, Vj) + Cj for j ^ i}, where ties are 
resolved in favor of the site with smaller additive weight. The regions V(vi, F + ) are called 
Voronoi cells. In the classic case, dist(-, •) has been defined as the Euclidean distance. Here, 
for dist(-, •), we use the weighted distance function c{y,x). 

Let v and v' be two points in F~ . We wish to study the intersection of the additive 
Voronoi diagram of v and v' with £ with respect to the weighted distance function. Without 
loss of generality, we assume that C = and C > 0, where C and C are the additive 
weights assigned to v and v', respectively. We denote the intersection of the Voronoi cell of 
v with £ by V(v, v', £; C), or simply by V(v) when no ambiguity arises. We have 

V(v,v',£;C) = V(v) = {x G £ : c(v, x) + C < c(v' , x)}. 

In Theorem 12. 11 we will show that if v and v' are at the same distance to F, then the Voronoi 
cell V(v) restricted to the line £ has a very nice structure (i.e., it is an interval). Furthermore, 
in Remark |2. 11 we will show that if v and v' are not at the same distance to F then Theorem 
l2.1l does not hold. In our algorithm, presented in Section 5, we use the information about the 
shape of V(v) in order to propagate approximate shortest paths in T> and it turns out that 
we need to only consider the sites that are restricted to be within a half-space of F and at 
the same distance to F . But, as we will see, this case in itself is mathematically challenging 
and provides valuable insights into the combinatorial structure of these diagrams. 

Theorem 2.1 The Voronoi cell V(v,v',£;C) is an interval on £ - possibly empty, finite or 
infinite. 

Proof: First, consider the case when C — 0. We denote the set of points x in F + such 
that c(i>,x) = c(i/,x) by B(v,v') and observe that it is a half-plane perpendicular to F. 
Therefore, the set of points x on t for which c(v,x) = c(v',x) is either a single point, the 
whole line I, or empty. Correspondingly, the Voronoi cell V(v,v',£;0) is either a half-line, 
empty, or the whole line £ and the theorem holds for C = 0. 

So, we assume that C > 0. We consider the equation c{v',x) — c(v,x) = C and claim 
that it cannot have more than two solutions. Before we prove that claim (Claim I2TT1 below) . 
we argue that it implies the theorem. 

Assume that the equation c(v', x) — c(v, x) = C has at most two solutions. If it does not 
have any or has just one solution, then the theorem follows straightforwardly. In the case 
where it has exactly two solutions, the cell V(v, v', £; C) has to be either a finite interval on I, 
or a complement of a finite interval on £. From the definition of the Voronoi cell V(v, v', £; C) 
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and C > it follows that V(v,v',£;C) C V(v,v' ,£;0). We know that V(v,v',£;0) is either 
empty, a half-line, or the whole line. If it is either empty or a half-line then V(v,v',£;C) 
must be either empty or a finite interval and the theorem holds. 

It remains to consider the case where V(v , v £; 0) is the whole line I. We argue that 
V(v,v',£;C) can not be a complement to a finite interval. We have V(v,v',£;0) = £ and 
therefore the line £ must be parallel to the half-plane B(v,v'). Furthermore, the plane 
containing B(v,v') is a perpendicular bisector of the segment (v, v ') and thus the point v' 
must have coordinates (0, y', z~) (Figure [2]). In this setting, using Lemma [2.1( e). we observe 
that c(v, x) and c(V, x) have same asymptotes at infinity and thus lim a ,_ 5 . 00 (c(t> / , x)—c(v, x)) = 

0. Therefore, the cell V(v,v',£;C) must be finite. The theorem follows. □ 
Next we establish the validity of the claim used in the proof of Theorem 12.11 

Claim 2.1 The equation c(v',x) — c(v,x) = C has at most two solutions. 

We will prove the claim by showing that the function g(x) = c(v',x) — c(v,x) is unimodal, 

1. e., it has at most one local extremum. We establish this property in two steps. First, we 
prove a characterization property of a local extremum of g (Proposition 12. ip . Then, we show 
that there may be no more than one point possessing that property. 

We focus our discussion on the case w~ ^ w + , since the other case is either simpler or 
can be treated analogously. We denote by a{x) and a'{x) the bending points defining the 
shortest paths from v and v' to x, respectively. We assume that £ is oriented and denote 
by et the positive direction unit vector on £. Furthermore, let a(x) and a'(x) be the angles 

between vectors v a(x), v'a'(x) and et, respectively (see Figure H]). 

These angles are completely defined by the angles if and 9 defining the corresponding 
shortest paths at the bending points a(x) and a'(x). Precisely, we have cos a = sin^cos^. 
Next, we prove that the angles a(xo) and a'(xo) must be equal at any local extremum xq. 

Proposition 2.1 If Xq is a local extremum of the function g, then a(xo) — at'(xo). 
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Figure 5: Illustration of the proof of inequality (J5J). 



Proof: The proof is by contradiction. Let us assume that a(x ) ^ a'(x ). We denote by 
(*k( x o) the angle between vectors a(x )x and e$. Similarly, a' K (x ) denotes the angle between 
vectors a'(xo)xo and et (Figure H]). The relation cos a = cos 9 sin (p and Snell's law readily 
imply that Kcosa K (xo) = cosa(xo) and kcosol' k (xo) = cosa'(xo), where k = w + /w~. Thus, 
we have that a K (xo) ^ a' K {xo). Without loss of generality, we assume that a K (xo) < a' K (xo). 
Under these assumptions, we show the existence of two points on xi and X2 on £, such that: 

x\ < xo < X2, g(xi) = g{x2)i and |a(x2)x 2 | + |a'(xi)xi| > |a(x2)xi| + |a / (xi)x 2 |. (4) 

By the assumption that xq is a local extremum of g, it follows that, for any positive real 
number 8 > inside the interval (xq — 8, x$ + 8), there are reals x\ and x^, such that 

xo — 8 < x\ < xq < x 2 < xq + 8 and g(x s x ) = g(x 2 ). 

On the other hand, if 8 converges to zero, then a(x 2 ) converges to a(xo), and of(xf) converges 
to a'(x ). Therefore, the inequality a K (x ) < o/ K (x ) implies that for a small enough 8 , the 
inequalities 

tt K (a ; i ) < a ' K { x5 \) an d Oi K {x 5 2 °) < c4( x 2°) 

hold. From these inequalities, it follows that if we make the quadrilateral a(x s 2 )a' {x{)x s 2 x{ 
planar by rotation of the point a(x 2 ) around £, then the obtained planar quadrilateral will 
be convex (Figure Therefore we have 

\a(x$)y$\ + |a'(zj°)xj°| > |a(4°) x i°l + W&i^l ( 5 ) 

which proves (jlj) for Xi = x^° and x 2 = x^ . 

Next, we estimate the sum c(v,xi) + c(v',X2) we show that there may be no more than 
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one point possessing that property. We use (J4]) and obtain 



c(v , xi) + c(v', x 2 ) = c(v, x 2 ) + c(V , Xi) 

= w~\va(x 2 )\ + u> + |a(x 2 )x 2 | + w~\v'a'(xi)\ + w + |a'(xi)x 1 | 
= w~(\va(x 2 )\ + \v'd(x-i)\) + w + (\a(x 2 )x. 2 \ + |a'(xi)xi|) 
> w~(\va(x 2 )\ + \v'a'(xx)\) + w + (|a(x 2 )x 1 | + |a'(xi)x 2 |) 
= w~\va(x 2 )\ + w + |a(x 2 )xi| + w~\v'a'(xi)\ + w + \a'(xi)x. 2 \. 

On the other hand, by the definition of the weighted distance function ([!]), we have the 
inequalities 



c(v, xi) < w \va(x 2 )\ + w + \a(x 2 )x.i\ and c(v',x 2 )<w \v'a'(xi)\ + w + \d(xi)x. 2 



that contradict the previous strict inequality Therefore, the angles a K (xo) and a' K (xo) and 
consequently a(xo) and a'(xo) must be equal. □ 
Our next step is to show that there cannot be two points on I satisfying Proposition 12.11 
To do that, we study in more detail the relationship between the position of points v and 
x, and the angle a(x). We observe that, for any fixed y (the ^-coordinate of v), there is a 
one-to-one correspondence between the real numbers x and the angles a. That is, for a fixed 
v, there is a one-to-one correspondence between the points x on £ and the angle between the 
shortest path 7f (v, x) and the positive direction on £. Hence, we may consider and study the 
well defined function x = x(y, a). We prove the following: 

Proposition 2.2 The second mixed derivative of the function x = x(y,a) exists and is 
negative, i.e., x ya < 0. 

Let us first show that Proposition 12.21 implies Claim 12.11 

Proof of Claim 1270 We assume that Proposition 12.21 holds and will show that the function 
g(x) has at most one local extremum. Recall that the point v has coordinates (0, y, z~) and 
let us denote the coordinates of the point v' by (x', y', z~). 

We first consider the case where y = y' . In this case, we observe that a'(x) = a(x + x'). 
In addition, the function a(x) is strictly monotone and, therefore, for any x, the angels a(x) 
and a'(x) are different. From Proposition 12.11 it follows that in this case the function g(x) 
has no local extremum. 

Next, we consider the case y ^ y'. We assume, for the sake of contradiction, that g(x) has 
two local extrema, say X\ and x 2 . By Proposition 12.11 a{x\) = a'(xi) and a(x 2 ) = a'(x 2 ). 
We denote ot\ = a{x\) = a'(xi) and a 2 = a(x 2 ) = a\x 2 ). Then, the difference x 2 — X\ can 
be represented, using the function x(y, a), in two ways 



x 2 - Xi 



ra2 pa-2 

I x a (y,a)da and x 2 — x\ = / x a (y',a)da. (6) 



Subtracting the last two equalities, we get 



OL2 ra2 ry poi2 

x a (y,a)da- / x a (y' ', a) da = / / x ya (y,a) dady. (7) 



y J ai 
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The integral on the right side is negative since by Proposition 12. 2\ the derivative x ya is 
negative, ot\ 7^ «2, and y 7^ y' . Hence, we have a contradiction and Claim |2~T1 follows. □ 

The proof of Proposition 12.21 is rather long and uses elaborate mathematical techniques 
and manipulations. On the other hand, the rest of the paper is independent of the details 
in that proof. So, we present the full proof for the interested readers in Appendix IA.1I 

Corollary 2.1 Consider a plane H in F + parallel to F and two points v and v' in F~ 
lying at the same distance from F. For any non-negative constant C , the Voronoi cell 
V(v, v', H; C) = {x G H : c(v, x) + C < c(v', x)} is convex. 

What can we say about the Voronoi cell of v in F + l The above corollary implies that the 
intersection between the Voronoi cell and any plane, parallel to F, is convex. This, as such, 
is not sufficient to conclude that the cell is convex. We close this section with the following 
conjecture. 

Conjecture 2.1 In the above setting, the Voronoi cell of v in F + , V(y , v', F + ; C) = {x G 

F + : c(v, x) + C < c(v', x)}, is convex. 

Remark 2.1 Examples showing that equal distance of the points v and v' from the bending 
plane F is a necessary condition for c(v', x) — c(v, x) to be unimodal in Theorem \2.1\ is not 
difficult to construct. In fact, if we take arbitrary points v and v' in F~ at different distances 
from F , it is very likely that c(v', x) — c(v, x) will have more than one local extrema and hence, 
for a proper choice ofC, the equation c(V, x) — c(v, x) = C will have more than two solutions. 
Such examples can be constructed by choosing arbitrary points v in F~ and x\, X2 on £ and 
then computing a point v' G F~ so that a(v,Xi) = a(v',Xi), for i = 1,2, where a's are the 
angles between the line i and the shortest paths coming from v and v' , respectively. As a 
result c(v',x) — c(v,x) will have local extrema at x\ and x%. 

3 Discretization of V 

In this section, we describe the definition of a carefully chosen set of additional points placed 
in V, called Steiner points. These Steiner points collectively form a discretization of T>, 
which is later used to approximate geodesic paths in T>. Steiner points are placed on the 
edges of T> and on the bisectors of the dihedral angles of the tetrahedra in T>. While it 
may seem more natural to place the Steiner points on the faces of the tetrahedra, placing 
them on the bisectors proves to be more efficient, leading to a speed up of approximately 
s~ x compared to the alternate placement. Recall that e is an approximation parameter in 
(0, 1). We provide a precise estimate on the number of Steiner points which depends on e 
and aspect ratios of tetrahedra of V. 

3.1 Placement of Steiner points 

We use the following definitions: 
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Definition 3.1 

(a) For a point x G T>, we define D(x) to be the union of the tetrahedra incident to x. We 
denote by dD(x) the set of faces on the boundary of D(x) that are not incident to x. 

(b) We define d(x) to be the minimum Euclidean distance from x to any point on dD(x). 

(c) For each vertex v G T>, we define a radius r(v) = d(v)/lA. 

(d) For any internal point x on an edge in T>, we define a radius r(x) = d(x)/24. The radius 
of an edge e EV is r(e) = max Ige r(x). 

Using radii r(v) and r(x) and our approximation parameter e, we define "small" regions 
around vertices and edges of T>, called vertex and edge vicinities, respectively. 

Definition 3.2 

(a) The convex hull of the intersection points of the ball B(v,er(v)) having center v and 
radius er(v) with the edges incident to v is called the vertex-vicinity of v and is denoted by 
D £ (v). 

(b) The convex hull of the intersections between the "spindle" \J xee B(x,er(x)) and the faces 
incident to e is called the edge- vicinity of e and is denoted by D e (e). 

On each edge e = AB of T>, we define a set of Steiner points as follows. Denote by AA 1 
and BB' the intersections of e with vertex vicinities D e (A) and D e (B), respectively. Points 
A' and B' are Steiner points. All other Steiner points on e are placed between A 1 and B'. 
Let M e be the point on e, such that d(M e ) = max iee cf(i). The point M e is defined to be a 
Steiner point. Next, we define a sequence of points M iy for i = 0, 1, . . . on M e A', by 



M = M e and |M;_iM;| = er(M) for i = 1,2 



All such points Mj between M e and A' are defined as Steiner points. Analogously, we define 
the set of Steiner points on M e B'. The number of Steiner points defined in this way on e is 
bounded by 

n 1 i 2 v. n 33 i \ AB \ 
C e -log-, where C e <- — log — ===== 

e e sma(e) ^Jr(A)r(B) 

and a(e) is the minimum angle between e and the faces on dD(e). (All logarithms with 
unspecified base are assumed to be base 2.) The total number of Steiner points placed on 
the edges of V is bounded by 61^ log |, where Ti is the average of the constants C e over 
all edges of V. 

The remaining Steiner points lie on the bisectors of the dihedral angles of tetrahedra in 
T>. Steiner points in any tetrahedron T of D are defined to lie on the six bisectors of the 
dihedral angles of T. Let the vertices of T be A, B, C, and D and let us consider one of 
the bisectors of the dihedral angles of T, say ABP (see Figure Efa)). Next, we describe the 
placement of Steiner points in the triangle ABP. 

Let the dihedral angle at AB of T be 7 and let PH be the height (altitude) of ABP (see 
Figure ISO 3 )). First, we define an infinite sequence of points Po, Pi, ... on PH by 

P = P, \Pi-iPi\ = y^j8\HPi\ sin( 7 /2), for % = 1, 2, . . . (9) 
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Then, we consider the sequence of lines Lj in the plane ABP, parallel to AB, and containing 
Pi, for % = 1,2,.... Let the intersection points of these lines with AP and BP be Ai and 
Bi, respectively. Points Ai and Bi lying outside of the vertex vicinities D £ (A) and D £ (B) 
are defined to be Steiner points, respectively. The intersection points of these lines with 
the boundary of the union of the edge vicinity D £ (AB) and vertex vicinities D £ (A) and 
D £ (B), are defined to be Steiner points. On each of the segments AiBi, we define a set of ki 
equidistantly placed Steiner points Pjj, j = 1, . . . , ki, where 

\AiBA 

and \PijP id+1 \ = 1 - , for j = 0,1,..., k t . (10) 

Ki + 1 

In the above expression, we have assumed that P^o = Ai and P^+i — Bi. 

Definition 3.3 The set of Steiner points in the triangle ABP consists of 

(a) all points Pi.j outside the union D £ (AB) U D £ (A) U D e (B), 

(b ) the intersection points of the lines Li with the boundary of that union. 

Next, we estimate the number of Steiner points placed in the triangle ABP. We denote 
h = Po = \PH\ and pi = \PiH\ for % = 1, 2, . . . In this notation, we have pi-\ — Pi = |PjPj_i| 
and |PjPj_i| = pi\Je/& sin (7/2), which implies 

p l = h\ i , where A = (1 + sin(7/2)) _1 . (11) 

Let i\ be the smallest index such that the line L i± is at distance smaller than er(e) from 
AB. We denote by K\ the number of Steiner points lying on lines Li, with i <i±, and by K 2 
the number of the remaining Steiner points in ABP. Let us estimate the number K\ first. 



\AiBi 



PiPi 



4-1 1 
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The number of Steiner points on a line Li, with i < ii, is ki + 2. Using (TTOT) and (fTTl) . we 
have 



A'; 







_ 1 PiPi+l 1 . 





(h- Pi )L4fl| 
- Pi+i) 



;i - \*)\ab\ 

hX^l - A) 



Thus, for the number we obtain 

ii-l 

5^(2 + fci) < 2(ii - 1) + 



^ 1 - A i 



i=i 



A' 



2(zi-l) + 



IABI 



1 - A* 1 . \ ^ nf . \AB\ 1 - A il_1 



(12) 



/i(l - A) \(1 — A)A il ^ 1 V - vx y /i(l-A) 2 A^- 1 
From the definition of ix and (jTTj) . we have hX n < er{e) < hX ll ~ l . Therefore, 



log. 



(13) 



(14) 



From ffl2l) and ffTB"]) . we obtain 

#1 < 7 yT 2 + 2 log A -i 

er(e)(l — A) 2 er{e) 

Next, we estimate K 2 , that is the number of Steiner points lying on segments AiBi with 
i>i\. By our definition, on the segment A^B^ there is a point M\, such that the triangle 
ABM\ lies entirely inside the edge vicinity. Let M 2 be the intersection point of the boundaries 
of the edge vicinity D £ (AB) and the vertex vicinity D £ (A) that lies in the triangle ABP. 
Similarly, let M' 2 ' be the intersection point of the boundaries of the edge vicinity D £ (AB), 
the vertex vicinity D e (B), and the triangle ABP. Furthermore, let i 2 be the smallest index 
such that the segment A^B^ is closer to AB than M 2 and similarly, let i 2 be the smallest 
index so that the segment AqB^ is closer to AB than M 2 . All Steiner points on segments 
AiBi, with i > ii, lie in the quadrilaterals A il A i ' 2 M 2 Mi and B^BqM^Mi. We denote the 
number of Steiner points in these two quadrilaterals by K 2 and K 2 , respectively. 

To estimate K' 2 , we show an upper bound on the number of Steiner points on AiBi, 
ii < i < i 2 that lie inside the quadrilateral A^A^ M 2 M\. Namely, if we denote this number 
by k[ and by Mj the intersection point between AiBi and AM\, we have, 

\AiMi\ \MxMx\pi , \A lx M x \ 



k[<2 



2 + 



Pi-Pi+x Ph(Pi-Pi+i) ' h\^(l-X)' 

Thus, the number of Steiner points inside the quadrilateral A^A^ M 2 Mi is bounded by 



K' 2 < 



(i' 2 - z0(2 + 



hX'i (1— A) 

quadrilateral Bi 1 BqM 2 Mi, we obtain K 2 < 
on K' 2 and K 2 , use (TTTT) . (fT3|) and obtain 



Analogously, for the number of Steiner points inside the 
— z'i)(2 + rTrTTT~TT *) • We sum the estimates 



(i' + i" 



2i x ) 2 + 



K 2 

\AB\a 



K' 2 + ^ < 
A* 1 ) 



H', + A. 



hX^{l — A) 



< 



& + if' 



2ii; 
-2ix) 



/lA^fl-A) 



2 + 



ZiA^(l-A) 

|AB| 
er(e)A(l - A) 



< 



(15) 
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From the definitions of the indices i' 2 , we easily derive that 

V2e 2 r(e)r(A) /.BAP V2e 2 r(e)r(B) ZPBA 

Pi' -l > r~TTTi cos , Pi" -i > r^m cos , 

FH \AM e \ 2 ' F * \BM e \ 2 ' 

where M e is the point on AB where the radius r(e) is achieved. These inequalities and (fTTj) 
imply 



■' ^ 1 _l i j •// / ! , i h\M e B\ 
t 9 < 1 + log A -i — = 75-777 and «o < 1 + log-v-i — = /nD . . 

v / 2£ 2 r(e)r(yl) cos 2 ~ A v / 2e 2 r(e)r(5) cos ZPBA 



Then, we use f fT3|) and obtain 



&IABI 



4 + «2 — 2i x < 2 + 2 log A -i - == === - 2 log A -i — 7^ = 

e 2 r(e) A /r(A)r( J B) er[e) 

2 + 2 log,-! 1 1 = 2 log,-, - 1 . 16 

e v /r(A)r( J B) e\y/r(A)r(B) 

Combining ( I14p . (|T5|) and (jlfip . we obtain 

IABI , IABI IABI , . 

+ 2 - 2 gr (e)A(l - A) Sa-1 e\^r(A)r(B) + er(e)(l - A) 2 

+ 21og A -i— p- + 41og A -i— 1 . 

From the last equation, we easily derive that 

K x + K 2 = C ABP (T)\ log-, 

where the constant Cabp{T) depends on the geometry of the tetrahedron T and is bounded 
by (see Appendix IA.2I for details) 

/ s \ AB \ , A\AB\ 2 h , , 

C ABP (T) < 23 1 , 1 log 1 18 
r(e)sm (7/2) r(e)r{A)r{B) 

Our discussion is summarized in the following lemma. 

Lemma 3.1 (a) The number of Steiner points placed on a bisector ABP of a dihe- 
dral angle 7 in a tetrahedron T, is bounded by Cabp(T)^ log - , where the constant 
Cabp{T) depends on the geometric features of T> around the edge AB and is bounded by 

\AB\ , 4\AB\ 2 h 



r(e) sin 2 ( 7 /2) U{ = r(e)r(A)r(B) ' 

(b) The number of segments that are parallel to AB on a bisector ABP, containing 
Steiner points, is bounded by C ABP (T)± log § , where C\ BP {T) < sin( ^ /2) log 2 r{ t)r{Ay{B) • 

(c) The total number of Steiner points is bounded by CiV) 1 ^ log where n is the number 
of tetrahedra in T> and C(T>) is the average of Cabp{T) over all 6n bisectors in T>. 
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By placing Steiner points in this way, in the next lemma, we show that it is possible to 
approximate the cell crossing segments that have their endpoints outside the vertex and the 
edge vicinities. 

Lemma 3.2 Let ABP be the bisector of a dihedral angle 7 formed by the faces ABC and 
ABD of a tetrahedron ABCD. Let X\ and x 2 be points on the faces ABC and ABD, 
respectively, that lie outside of the union D £ (AB) U D e (A) U D £ (B). Then, there exists 
a Steiner point q on ABP, such that max(Zx 2 x\q, Zx\x 2 q) < \J\ and \xiq\ + \qx 2 \ < 
(l + e/2)\ Xl x 2 \. 

Proof: Clearly, the segment Xix 2 intersects the bisector triangle ABP in a point Xq lying 
outside the vertex vicinities D e (A), D E (B), and the edge vicinity D e (AB). Recall that 
Steiner points in ABP are placed on a set of lines Lj parallel to AB and passing through 
the sequence of points Pi on the altitude PH of ABP. Let i be the maximum index such 
that the line L io is farther away from AB than from xq. We define q to be the closest Steiner 
point to xq on the line Lj . 



D 




B 



Figure 7: Illustrates Lemma [3.21 

First, we estimate the angles Zx 2 x±q = Zx Xiq and /Lx\X 2 q = ZxoX 2 q- By our definition 
of the Steiner points and Pythagorean theorem, it follows that 

|x g|<^A^y|sin|, (19) 

where h and A are as defined above (see pip ). Let p be the radius of the smallest sphere 
containing xq and q and touching the face ABC. It is easily observed that 

2p > (hX k) - M) sin 1 > 8v ^-/^ A* sin 2 . (20) 
' V 2^28^ 2 v 1 
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If we denote the angle /Lx§X\q by 9\, then sin^i < ^7^, and using (fT9|) and ([20]) . we obtain 

sin ()l <M < y|. (21) 

The same estimate applies to angle Z-X^x^q- Hence, the first inequality of the lemma holds. 
Next, we prove the second inequality. We denote by 9, 9±, and 9% the angles of the triangle 
qx\X2 at q, x\ and X2, respectively (Figure Ej). By a trigonometric equality valid in any 
triangle, we have 

V sin(0/2) J 

Thus, it suffices to prove that 2&m ( 6l i^W( 6 '2/ 2 ) < e ji. By (|2"Tj) . it follows that sin^x and sin# 2 
are smaller than \fsj2 and from e < 1 we have 6* > 7r/2. Therefore, we obtain 

2 sin(6»i/2) sin(# 2 /2) sin 9 X sin # 2 £ 



sin(0/2) 2 sin(0/2) cos(^/2) cos(fl 2 /2) ~ 4 sin(0/2) cos(0i/2) cos(# 2 /2) 



e ee 

< < -•: : 



4sin(0/2)(sin(0/2) +sin(0 1 /2)sin(0 2 /2)) ~ 4sin 2 (0/2) ~ 2' 



4 Discrete paths 

In this section, we use the Steiner points introduced above for the construction of a weighted 
graph G £ = (V(G £ ), E(G £ )). We estimate the number of its nodes and edges and then 
establish that shortest paths in T> can be approximated by paths in G £ . We follow the 
approach laid out in [4], but the details are substantially different, as we have to handle 
both the vertex and edge vicinities, as well as the bisectors in 3-d space. 

The set of nodes V(G £ ) consists of the vertices of V, the Steiner points placed on the 
edges of V and the Steiner points placed on the bisectors. The edges of the graph G £ join 
nodes lying on neighboring bisectors as defined below. A bisector is a neighbor to itself. Two 
different bisectors are neighbors if the dihedral angles they split share a common face. We 
say that a pair of bisectors sharing a face / are neighbors with respect to /. (So, a single 
bisector b is a neighbor to itself with respect to both faces forming the dihedral angle it 
splits.) 

First, we define edges joining pairs of Steiner points on neighboring bisectors. Let p 
and q be nodes corresponding to Steiner points lying on neighboring bisectors b and bi, 
respectively, that share a common face /. We consider the shortest weighted path between 
p and q of the type {p, x, y, q}, where x and y belong to / (points x and y are not necessarily 
different). We refer to this shortest path as a local shortest path between p and q crossing / 
and denote it by Tt{p, q; /). Nodes p and q are joined by an edge in G £ if none of the points 
x or y are on an edge of /. Such an edge is said to cross the face /. In the case where p and 
q lie on the same bisector, say b, splitting an angle between faces f\ and / 2 , we define two 
parallel edges in G £ joining p and q - one crossing f\ and another crossing / 2 . 
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The cost of an edge (p, q) in G £ that crosses a face / is defined as the cost of the local 
shortest path 7r(p, q; f) and is denoted by c(p, q; f), or simply by c(p, q) when no ambiguity 
arises. Formally, we have 

c(p,q) = c(p,q;f) = \\it(p,q; f)\\ = min(||px|| + \\xy\\ + \\yq\\). (22) 

x,yef 

Next, we consider a node p of G £ lying on an edge e of T>. The node p can be either a 
Steiner point on e or a vertex of T> incident to e. It is adjacent to nodes lying in tetrahedra 
in D(e). The edges of G £ incident to p are associated with pairs of neighboring bisectors as 
follows. We consider a tetrahedron t in D(e), and describe edges incident to p in t. Let fx 
and f2 be the two faces of t incident to e, and let b be the bisector of the dihedral angle 
formed by j\ and fi- We define edges between p and nodes lying on bisectors in t that are 
neighbors of b. There are four such bisectors - two with respect to f\ and two with respect 
to f 2 . For a node q on a neighboring bisector b 1 sharing, say, the face fi with b, we consider 
the local shortest path 7r(p, q;fi). By definition, Tt(p,q;fi) = {p,x,q}, where x G f\. We 
define an edge between p and q if and only if the point x defining the local shortest path is 
in the interior of f\. The cost of the edge (p, q) equals the cost of the local shortest path 
7r(p,9;/i), i-e-, 

c(p,q) = c(p,q;fi) = \\ir(p,q;fi)\\ = mm(\\px\\ + \\xq\\). 

We associate the edge (p, q) to b, bx and j\ and say that it crosses f\. Furthermore, p is 
joined to nodes on b by pair of parallel edges, provided that the corresponding local shortest 
paths do not touch the edges of T> - one crossing fi and the other crossing / 2 . 

Lemma 4.1 We have \V{G e )\ = O(^logi) and \E(G £ )\ = 0(^log 2 i). 

Proof: The estimate on the number of nodes follows directly from Lemma [3.11 and the fact 
that T> has 0(n) vertices. The number of edges in G e can be estimated as follows. There 
are 0{n) faces in T> and at most 21 pairs of neighbor bisectors with respect to a fixed face 
in T>. By Lemma l3TTl (a). there are 0{\ log 2 -) pairs of nodes lying on two fixed neighboring 
bisectors. When combined, these three facts prove the estimate on the number of edges of 
G £ . □ 
Paths in G £ are called discrete paths. The cost, c(ir), of a discrete path 7r is the sum of 
the costs of its edges. Note that if we replace each of the edges in a discrete path 7r by the 
corresponding (at most three) segments forming the shortest path used to compute its cost 
we obtain a path in T> with cost c(ir). Next, we state the main theorem of this section. 

Theorem 4.1 Let tc(vo,v) be a shortest path between two different vertices vq and v in T>. 
There exists a discrete path n(v ,v), such that c(ii(vo,v)) < (1 + e)||ff(i;o, 

Proof: We prove the theorem by constructing a discrete path tt(vo,v) whose cost is as 
required. Recall that the shortest path tt(vo,v) is a linear path consisting of cell-crossing, 
face-using, and edge- using segments that satisfy Snell's law at each bending point. We 
construct the discrete path ir by successive modifications of n described below as steps. 
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(a) 



(b) 



Figure 8: (a) Replacement of a cell-crossing segment s = (21,22) by a two-segment path 
{xi,p, 22}. (b) Replacement of fti by a path n 2 joining Steiner points p^. Note that edges 
(PiPi+i) denote local shortest paths, rather than straight-line segments. 

Step 1: In this step, we replace each of the cell-crossing segments of tt, which satisfy the 
conditions of Lemma 13.21 by a two-segment path through a Steiner point. Precisely, let 
s = (21,22) be a cell-crossing segment in n (Figure M (a))- Let fi and / 2 be the faces 
containing 21 and 22, respectively. Let e = (A, B) be the common edge between fi and 
/ 2 . Assume that s is outside of the union of the edge and vertex vicinities D e (e) U D e (A) U 
D £ (B). We refer to such segment as vicinity-fre^. Then, according to Lemma r3.2[ there is a 
Steiner point p on the bisector b splitting the dihedral angle formed by j\ and / 2 such that 
\x±p\ + \px 2 \ < (1 + e/2)\xix 2 \. So, in this step, each cell-crossing and vicinity-free segment 
s = (21, 2 2 ) is replaced by two-segment path {x\,p, 2 2 }, where p is the approximating Steiner 
point as described above. Clearly, after this step, we obtain a path joining vo and v, whose 
cost does not exceed (1 + e/2)||7f||. We denote this path by fri (see Figure [8] (b)). 

Step 2: In this step, we consider the sequence of Steiner points added as new bending points 
along 7Ti in Step 1. In the case where two consecutive Steiner points are split by a single 
bending point or a face-using segment on tt\, we replace the sub-path between them by the 
corresponding local shortest path. Precisely, assume that p 1 and p 2 are consecutive Steiner 
points along fri and the sub-path between them is either {p±,x,p 2 } or {pi,x,y,p 2 }, 2 and 
y are bending points on the face /, shared by the two neighboring tetrahedra containing p\ 
and p 2 , respectively. So, in Step 2, we replace all such sub-paths by the local shortest paths 
7r(Pi)P2; /) = {pi,x,y,p 2 }, using (122]) . We denote the obtained path by tt 2 (Figure [8] (b)). 
Clearly, n 2 is a path joining v and v, whose cost does not exceed that of 7Ti. Hence, 



In the following two steps, we identify the portions of tt 2 that lie inside the vertex and 
edge vicinities and replace them with discrete paths using the corresponding vertices and 
edges. 

4 Note that such a segment still can have an end-point in a vertex or edge- vicinity related to other vertices 
or edges incident to fi and fa- 



7T 2 || < INI < (l + e/2)||7f 



(23) 
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Step 3: Follow tt 2 from vq to v and let be the last bending point on tt 2 that lies inside the 
vertex vicinity D E (v ). Next, let bi be the first bending point after a that is in the vertex 
vicinity, say D e {vi). Likewise, let a± be the last bending point in D £ (vi). Continuing in this 
way, we define a sequence of, say k + 1 for some k > 1, different vertices Vq, V\, . . . , Vk — v and 
a sequence of bending points do, 6i, ai, . . . , a^-i, on 7r 2 , such that for % — 0, . . . , k, points 
bi,cii are in D £ {vj) (we assume b$ = Vq, = v). Furthermore, by our definition, portions of 
7r 2 between and b i+1 do not intersect any vertex vicinities. We partition n 2 into portions 

^2^0, do), 7f 2 (a , &l), 7T 2 (6i, «1 ) , • • • , 7T 2 (6 fe , «)• (24) 

The portions 7r 2 (aj, 6j+i), for i = 0, . . . , k — 1, are called the between-vertex-vicinitiesportions, 
while the portions 7f 2 (6i, a*), for z = 0, . . . , fc, are called the vertex-vicinity portions. 

We define path 7f 3 by replacing each of the vertex-vicinities portions by a two segment 
path trough the corresponding vertex and show that the cost of 7r 3 is bounded by (1 + 
e - / 6) 1 1 -7T2 1 1 . Consider a between-vertex-vicinities portion 7r2(ai,&i+i) for some < i < k — 1. 
If this portion consists of a single segment (aj, &t+i), then the vertices Vi and must be 
adjacent in V and we define 7r 3 (wj, Ui+i) to be the segment (utjUi+i). The length of (vi,v i+ i) 
is estimated by using the triangle inequality and the definition of the vertex-vicinities as 
follows: 

\viV i+1 \ < \vidil + |ai& i+ i| + \b i+1 v i+1 \ < \aib i+1 \ + e(r(vi) + r(v i+1 )) < 

\aA+i\ + —( d ( l, i) + d i v i+l)) < + ~\ v i v i+l\- ( 25 ) 

To estimate the cost of the segment (v v i+ i), we observe that (a i5 lies inside a tetrahe- 
dron incident to (vi,v i+ i). Thus, the weight of (vi,v i+ i) is at most the weight of (ai,b i+ i). 
This observation and ( |25l) readily imply 

11^3(^,^+1)11 = ||viVi+i|| < (l + |)||ai6 i+1 || = (1 + |)||7r 2 (ai,6i+i)||- (26) 

o o 

In the general case, where ^(oj, b i+ i) contains at least two segments, we follow the bending 
points along 7r 2 (ai, b i+1 ) and define X to be the last bending point on the boundary dD(yi) 
(see Definition 13. 1 f) . If the path 7r2(aj,6j + i) lies entirely in D(vi), then we set X = 644.1. 
Thus, the bending points on the path 7T2 between a, and X lie in the tetrahedra incident 
to V{. Let to, be the minimum weight among the segments of the path 7i2(ai,X) and let x 
be the first bending point after a, incident to a segment, whose weight is Wi. Analogously, 
define the bending points Y and y, by following the bending points of the backward path 
7r 2 (6j4i,aj) from 6 i+1 . Note that x precedes y on the path 7f 2 (ai, 6 i+1 ). We define the path 
fr 3 (vi,Vi + i) as the concatenation of the segments (vi,x), (y,v i+ i) and the portion 7r 2 (x,y), 
i.e., 

n 3 (vi, v i+1 ) = {(v u x), tt 2 (x, y), (y, v i+1 )}. 

Next, we estimate the cost of 773(^,^4.1). First, we observe that the weight of the segment 
(vi,x) cannot exceed u>j. Then, we use the triangle inequality and the fact that ai is inside 
the vertex vicinity D e (vi), obtaining 

\\vix\\ < Wi\viai\ + Wi\n 2 {ai,x)\ < Wi\viai\ + \\TT 2 {ai,x)\\ < 
Wier{vi) + ||7r 2 (aj,x)|| < Wi^-d{vi) + \\ii 2 (a u x)\\. 
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Analogously, for the cost of the segment (y,v i+ i), we have 



yv i+ i\\ < w i+1 —d{v i+ i) + ||7f 2 (j/,6i+i)||- 



Using these estimates, and the way we defined the path ^(ui, i>i+i), the weights Wi, Wi+i, 
the distances d(vi), d(vi + i), and the points X, Y, we obtain 

||7r 3 0i,w i+ i|| < ||7r 2 (ai,6i+i)|| + -^(widfa) + w i+1 d(v i+1 )) < (27) 

£ £ 

((^(ai^i+i)!! + — (||vr 3 (^,X)|| + ||7f 3 (F,^ + i)||) < ||7r 2 (ai,6i+i)|| + j(\\^3(vi,v i+ i) ||, 

which implies the estimate (1261) in the general case. Applying the above construction to each 
pair of consecutive vertices in the sequence v , Vi, . . . , v k = v , we obtain a linear path 

^3(^0, V) = {n 3 (v , Vi), 7T 3 (fl, V 2 ),- ■ .,TT 3 {v k -l,v)}, 

that has no bending points inside vertex vicinities except for the vertices vq, V\, . . . , v k = v . 
We estimate the cost of this path by summing up fT2"6"j) . for i — 0, . . . , k — 1, and obtain 



Observe that the path 7r 3 constructed above may contain self intersections (e.g., if one and 
the same vertex vicinity is visited twice by 7r 2 ) . It is also possible that tt 3 may contain 
consecutive face-using segments. Hence, at the end of Step 3, we traverse the obtained path 
and compress it. That is, we remove the loops in case of self intersections. We replace the 
consecutive face-using segments (which obviously lie in the same face) by the single face- 
using segment joining their free end-points. We denote the compressed path again by 7f 3 . 
Clearly, compressing reduces the cost of the path and hence the estimate ( )28|) remains true 

for 7T 3 . 

Next, in Step 4, using a similar approach as above, we further partition each vertex- 
vicinity-portion vr 3 (i>i, fj+i) into between- edge-vicinities portions and edge-vicinity portions. 
Then, we replace each edge-vicinity portion by an edge-using segment plus 2 additional 
segments and estimate the cost of the resulting path 7r 4 . 

Step 4: First we define analogues of vertex and between-vertex vicinities for edges. Let 
(vi,a) be the first segment of the path ^3(^,^+1). If a is not inside an edge-vicinity, 
we define a^o = v%. Otherwise, if a is inside an edge- vicinity, say D £ (eo), and let a' 
be the first bending point on the path ^3(1^,1^+1) after Vi lying on dD(eo), then we de- 
fine a^o to be the last bending point on n 3 (vi,a') that is inside D e (eo). Next, let b^i 
be the first bending point on 7r 3 (aj j0 , v j+i) that is inside an edge-vicinity, say D e (e\) and 
let b' be the first bending point on ^3(61,1,^+1) that is on dD(ei). We define a^i as 
the last bending point on 7r 3 (6ji,6') that is in the same edge vicinity as b^i. Assume 
that, following this approach, the sequence of bending points a i>0 , b it i, a^i, . . . , a i>ki -i, b itki 
has been defined. They partition the portion vr 3 (fi,fj+i) into sub-portions n 3 (vi,v i+ i) = 



k-l 




(28) 
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Figure 9: Replacement of subpath 7ts(bij, a it j) by the edge-using segment (qijPij). 



{# 3 (t>j, a ii0 ), ... ,^3(0^-1 Portions between a iyj and 

for j = 0, . . . , ki — 1, are called the between-edge-vicinity portions. Portions between 
b it j and a it j, for j = 0, . . . , fcj, are called the edge-vicinity portions (b i = Vi and = v i+ i). 

According to our construction, the bending points a^o, b^x, a^±, . . . , a^-i, 6^, defining 
the above partition lie inside edge vicinities. Moreover, consecutive points fejj and Ojj, for 
j = 0, . . . , ki, are in one and the same edge-vicinity D E (e.j). 

For j = 0, . . . , let b\ • and be the orthogonal projections of the points fty and 
dij onto the edge e^, respectively (Figure E]). Let p i} j and q^ be the Steiner points on ej 
defining the largest sub-interval of the interval (a^ ■, 6^) on and assume that p^j is between 
a^ - and q^j. (In the case where the interval (a' ij -,6^ J ) contains no Steiner points, we define 
Pi,j = %j to be the closest Steiner point to on ej.) In 7r 4 , the edge-using segment (qijPij) 
will replace in 7r 3 the subpath 7^3(6, j, Let us estimate the resulting error. From the 

definition of the edge vicinity D E (ej), the Steiner points on e^, and the radii r(p it j) and r{q^j), 
it is easy to derive that 

3 3 

\Pi,jdi,j\ < ^r(pij) and l^^jl < -er(q i;j ). (29) 

Furthermore, by our construction and the fact that (qij,Pi,j) is an edge-using segment, it 
follows that 

llftjPijII < \\^3(bi,j,aij)\\. (30) 
Next, we modify between-edge- vicinities portions 713(0^, fcij+i), into paths 7f 4 (pjj, ftj+i), 
joining Steiner points pjj and qij+i, and not intersecting any vertex or edge vicinities. We 
apply a construction analogous to the one used in Step 3 to define the paths n 3 (vi,v i+ i). 

We fix j and consider the between-edge- vicinities portion 773(0^, bij+i). We first consider 
the special case where 713(0^-, frij+i) is the segment (a^j, feij+i). In this case, we observe that 
ej and e J+ i must be edges of the tetrahedron containing the segment (ajj, fcij+i) and define 
n^pij, qij+i) = (pi,j,Qi,j+i)- We estimate the length of this segment using the triangle 
inequality, the estimates ff29l) and Definition 13.11 as follows 

\Pi,jqi,j+l\ < \Pi,j a i,j\ + \ a i,jbi,3+l\ + < y( r fe,j) + r (<?M+l)) + KA.j+l| < 

^( d (P<j) + d (<li,j+i)) + KAi+i| < ^\Pijqi,j+i\ + \aijb iJ+1 \. 

Using this estimate and the observation that the weight of the segment (pij,Qij+i) cannot 
exceed the weight of (a^j, we obtain 

\\MPi,j,<li,3+l)\\ = \\Pi,3<lij+l\\ < (1 + |)lk,Aj+ill = (1 + ^II^K^Aj+i)!!- (31) 
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Next, we consider the general case, where ||7r 3 (ajj, consists of at least two seg- 

ments. Let X be the first bending point after a it j that is on the boundary dD(ej). If the 
path 7f 3 (ajj, fejj+i) is entirely inside D(ej), then we set X = Furthermore, let w'j be 

the minimum weight among the segments in 7r 3 (ajj, X), and let x be the first bending point 
after that is an end-point of a segment whose weight is w' y We define the weight it^- +1 , 
and the bending points Y, and y, analogously with respect to ftij+i and the edge vicinity 
D e (ej + i). It follows that the point x precedes y along 7T 3 (ajj, We define the portion 

of the path vr 4 joining p itj and g ij+ i by 7f 4 (Pij, ftj+i) = {(Pij, a;), 7r 3 (:r, 2/), (?/, and 
estimate its cost. Let us first estimate the cost of the segment (pij,x). We observe that 
Ibij^ll < ^bij x l) that \pij,x\ < \pijCLij\ + |#3(ai,j, x)\ (by triangle inequality), and that 
the segments on the path n 3 (aij,x) have weight greater than or equal to w' y Using these 
observations, (j29p . and Definition 13. 1\ we obtain 

\\Pijx\\ < w'j\Pij<*ij\ +^|7f 3 (a ijj ,z)| < Wj\pi d a i;j \ + Wn^Oij^x^ < 



3e . , 



™'j r (Pi,j) + \\^3{ai,j,x)\\ = —w'jdipij) + \\TT 3 (a id ,x)\\. (32) 



Analogously, we have 



\yqi,j+i\\ < ^' j+1 d{q iJ+1 ) + IMl/Aj+OH- (33) 



Using the definition of the path 7r 4 (pij, qij+i), the estimates ( 1321) and (1331) . and the definition 
of the distances d(pij) and <i(g ij+ i), the weights w'j and w'j +1 , and the points X and Y, we 
obtain 

IMPij, qi,j+i)\\ = \\M a i,3, kj+i)\\ + Y^(^j d (Pij) + v>' j+1 d(q iJ+1 )) (34) 

< \\f3{%j,b i)j+1 )\\ + ^(ii^^x)!! + ||7f 3 (y,^- +1 )||) 

< ||7f 3 K,i,&i,i+i)|| + ^\\^3(pi,j,qi,j+i)\\, 

which implies estimate (13 lj) . in the general case. Finally, combining segments (qi,j,Pi,j) and 
paths fr^pij, qij+i), for j = 0, . . . , /cj, we construct a path 

^4(^,^+1) = {(ui,Pi,o), ^4(^,0,91,1), (%,i,Pi,i), ^4^,1,^,2), • • • ,7r4(Pi,fe i -i,9i,feJ, (ft,**, Ui+i)}- 

This path has no bending points in any of the edge or vertex vicinities. Its cost can be 
bounded using ( 130]) and (j3lj) as follows 

fc j 1 

||7f 4 (^,U i+ i)|| = ^2\\qi,j,Pi,j\\ + ^2\\TT4(Pi,j,qi,j+l)\\ < (35) 

i=o j=o 

hi hi 1 



\\M b ij7<kj)\\ + 5^(1 + ^II^Kj'Aj+OII < (1 + j)\\^3(Vi,V i+1 ) 



e 

U + 

i=o 3=0 
where we assume t»j = b i>0 = q ifi and v i+1 = a iikl = Pi,k 
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The paths #4(1^, v i+ i), for i = 0, ...k — 1, form a linear path 7r 4 (w ,?;), whose cost is 
estimated using ( I3"5j) . (|2"B"j) . and ( 12"B"|) by 

fe-i 

7f4(V0,«) < ^(l + |)||^3(^,^ + l)ll = (l+y)|Mvo,T>)|| < (1+|)II^(^0^)|| < (l+e)||7f(V0,v)||. 
i=0 

(36) 

As in Step 3, it is possible for 7r 4 to contain self-intersections and consecutive face-using 
segments. Hence, we traverse 7r 4 and compress it by removing loops and by replacing con- 
secutive face- using segments. The obtained path is denoted again by 7r 4 , and estimate (|36|) 
is valid. 

The bending points defining 7r 4 can be partitioned into two groups. The first group 
consists of bending points corresponding to nodes of the graph G £ , i.e., Steiner points on 
bisectors, Steiner points on edges, and vertices of T>. The second group consists of the 
remaining bending points of 7f 4 , which are bending points inside the faces of T>. We complete 
the proof of the theorem by showing that the sequence of the nodes in the first group defines 
a discrete path tt(vo,v) whose cost c(7r(t> , v)) < \\tt4(vo,v)\\. It suffices to show that any two 
consecutive nodes (bending points in the first group) along the path 7r 4 are adjacent in the 
approximation graph G £ . 

To show this, we review closely the structure of the path 7r 4 . In Step 3, portions of #2 
related to vertex vicinities have been replaced by two segment portions through-vertices of 
T>. Furthermore, we observe that the segments (vi,x) created in Step 3 are either a face- 
using segments or join Vi to a Steiner point on a bisector. The same applies to segments 
(y, Vi + i). Similarly, in Step 4, portions related to edge- vicinities have been replaced by three 
segment portions visiting corresponding edges. Again segments (pij,x) are either face-using 
segments or join p it j to a node on a bisector, that is a neighbor of the bisector incident to the 
edge containing p^j. The same applies to the segments (y, q^j+i)- In summary, the segments 
created in Steps 3 and 4 are of one of the following two types: 

1. A face-using segment with one of its endpoint being a (node) vertex of T> or a Steiner 
point on an edge of T>. 

2. A segment joining two nodes, at least one of them being a Steiner point on an edge of 
Dora vertex of T>. 

The remaining segments in 7r 4 are cell-crossing and face-using segments, whose endpoints are 
outside any vertex or edge vicinity. All the cell-crossing segments in 7r 4 were created during 
Steps 1 and 2. Hence, one of their endpoints is a (node) Steiner point on a bisector of a 
tetrahedron. Finally, due to the compressing, there are no consecutive face-using segments 
in 7f 4 . 

Now, let p and q be two consecutive nodes along the path 7r 4 . We show that p and q 
are adjacent in G £ . We consider, first, the case where at least one of the nodes, say p, is a 
vertex of T>. Let x be the bending point following p along the path 7r 4 . By the definition 
of bending points adjacent to the vertices (in Step 3), we know that (p, x) is a face-using 
segment followed by a cell-crossing segment (x, xi), joining x to a (node) Steiner point on a 
bisector lying in one of the tetrahedra incident to the face that contains (p, x). So, q = X\ 
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and q is inside a tetrahedron incident to p. Thus, p and q are adjacent in G £ . The case where 
at least one of the nodes p or q is a Steiner point on an edge of V can be treated analogously. 

Assume now that both p and q are Steiner points on bisectors. Let x and X\ be the 
bending points following p along 714. The point x has to be a bending point on a face of 
the tetrahedron containing p. The segment (x, x\) is either a cell-crossing or a face-using 
segment. In the first case, q must coincide with x\ and is adjacent to p in G £ , since it lies in 
a tetrahedron that is a neighbor to the one containing p. In the second case, where (x,Xi) 
is a face- using segment, we consider the bending point x 2 that follows Xi along the path 7r 4 . 
The segment (xi,x 2 ) must be a cell-crossing segment. Thus, in this case, q = x 2 is adjacent 
to p, because the tetrahedra containing p and q are neighbors. 

We have shown that any pair p and q of consecutive nodes on the path 714 are adjacent in 
G £ . Hence, we define a discrete path n(v , v) to be the path in G e following the sequence of 
nodes along 7r 4 . Finally, we observe that the sub-paths of n^p, q) joining pairs of consecutive 
nodes stay in the union of the tetrahedra containing these nodes and cross faces shared by 
the bisectors containing them. Hence, by the definition of the cost of the edges in G £ , we 
have c(p,q) < I^Qo, q)\\. Summing these estimates, for all edges of ti(vq,v), and using ([36]). 
we obtain, c(ir(v , v)) < \\tt 4 (v ,v)\\ < (1 + e)\\n(v ,v)\\. □ 

5 An algorithm for computing SSSP in G £ 

In this section we present our algorithm for solving the Single Source Shortest Paths (SSSP) 
problem in the approximation graph G e = (V(G e ), E(G £ )). Straightforwardly, one can apply 
Dijkstra's algorithm, which runs in 0(\E(G £ )\ + |V((j e )| log ^(G^)!) time. By Lemma [4.11 
we have \V{G e )\ = 0(^log±) and \E(G £ )\ = 0(fJog 2 \). Thus, the SSSP problem in G £ 
can be solved in 0{\ log - log -) time. 

In the remainder of this section, we demonstrate how geometric properties of our model 
can be used to obtain a more efficient algorithm for solving the SSSP problem. More precisely, 
we present an algorithm that runs in O ( | | (log \ V £ \ + log 3 ^)) = 0{-^ log j log 3 time. 

Informally, the idea is to avoid consideration of large portions of the edges of the graph G £ 
when searching for shortest paths. We achieve that by applying the strategy proposed first 
in [25l [26] and developed further in [I] and by using the properties of the weighted distance 
function and additive Voronoi diagrams studied in Section I2~2l We maintain a priority queue 
containing candidate shortest paths. At each iteration of the algorithm, a shortest path from 
the source s to some node u of G £ is found. Then, the algorithm constructs edges adjacent to 
u that can be continuations of the shortest path from s to u and inserts them in the priority 
queue as new candidate shortest paths. In general, one needs to consider all edges adjacent 
to u as possible continuations. In our case, we divide the edges adjacent to u into log -) 
groups related to the segments containing Steiner points in the neighboring bisectors and 
demonstrate that we can consider just a constant number of edges in each group. The latter 
is possible due to the structure of the Voronoi cell V(u) of the node u in the additive Voronoi 
diagram related to a fixed group (see Theorem 12. ip . 

This section is organized as follows: In the next subsection, we describe the general 
structure of the algorithm. In Subsection 15.21 we show how this strategy can be applied in 
our case and present an outline of the algorithm. We provide details of the implementation 
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Figure 10: The sets S, S a , and an optimal edge e(S) = (u,v) in E(S) are illustrated. A shortest 
path from s to u is illustrated by a dashed curve. 

of the algorithm and analyze its complexity. Finally, at the end we establish the main result 
of the paper. 

5.1 General structure of the algorithm 

Let G(V, E) be a directed graph with positive costs (lengths) assigned to its edges and s 
be a fixed node of G, called the source. A standard greedy approach for solving the SSSP 
problem works as follows: a subset, S, of nodes to which the shortest path has already been 
found and a set, E(S), of edges connecting S with S a C V \ S are maintained. The set S a 
consists of nodes not in S but adjacent to S. In each iteration, an optimal edge e(S) = (u,v) 
in E(S) is selected, with source u in S and target v in S a (see Figure [TO]) . The target vertex 
v is added to S and E(S) is updated correspondingly. An edge e = e(S) is optimal if it 
minimizes the value 5(u) + c(e), where S(u) is the distance from s to m and c(e) is the cost 
of e. 

Different strategies for maintaining information about E(S) and finding an optimal edge 
e(S) during each iteration result in different algorithms for computing SSSP. For example, 
Dijkstra's algorithm maintains only a subset Q(S) of E(S), which, however, always contains 
an optimal edge. Alternatively, as in [3], one may maintain a subset of E(S) containing 
one edge per node u in S. The target node of this edge is called the representative of u 
and is denoted by p(u). The node u itself is called predecessor of its representative. The 
representative p(u) is defined to be the target of the minimum cost edge in the propagation 
set I{u) of u, where I{u) C E(S) consists of all edges (u,v) such that 5(u) + c(u,v) < 
S(u') + c(u',v) for all nodes u' G S that have entered S before u. The union of propagation 
sets forms a subset Q(S) of E(S) that always contains an optimal edge. Propagation sets 
I(u), for u G S, form a partition of Q(S). The propagation sets of the vertices in S form a 
partition of E(S), which is called propagation diagram, and is denoted by Z(S). 

The set of representatives R C S a can be organized in a priority queue, where the key of 
the node p{u) in R is defined to be 5{u) + c(u,p(u)). Observe that the edge corresponding 
to the minimum in R is an optimal edge for S. In each iteration, the minimum key node v 
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in R is selected and the following three steps are carried: 

Step 1. The node v is moved from R into S . Then, the propagation set I(v) is computed 
and the propagation diagram T(S) is updated accordingly. 

Step 2. Representative p(v) of v and a new representative, p(u), for the predecessor u of v 
are computed. 

Step 3. The new representatives, p{u) and p{v), are either inserted into R together with 
their corresponding keys, or ( if they are already in R) their keys are updated. 

Clearly, this leads to a correct algorithm for solving the SSSP problem in G. The total 
time for the priority queue operations is 0(|V| log |V|). Therefore, the efficiency of this 
strategy depends on the maintenance of the propagation diagram, the complexity of the 
propagation sets, and the efficient updates of the new representatives. In the next subsection, 
we address these issues and provide necessary details. 

5.2 Implementation details and analysis 
5.2.1 Notation and algorithm outline 

Our algorithm follows the general strategy as described in the previous subsection. First, 
we convert G e into a directed graph by replacing each of its edges by a pair of oppositely 
oriented edges with cost equal to the cost of the original edge. 

Let, as above, S be the set of the nodes to which the shortest path has already been 
found and E(S) be the set of the edges joining S with S a C V\S. We partition the edges of 
G e (and respectively E(S)) into groups so that the propagation sets and the corresponding 
propagation diagrams, when restricted to a fixed group, have a simple structure and can be 
updated efficiently Then, for each node u in S, we will keep multiple representatives in R - 
a constant number on the average, for each group where edges incident to u participate and 
where its propagation set is non-empty. A node in S a will have multiple predecessors - at 
most as many as the number of the groups where edges incident to it participate. We will 
show that the number of the groups, where edges incident to u can participate, is bounded 
by O(-^log^) times the number of bisectors incident to u. In a fixed group, we will be 

able to compute new representatives in O(log-) time and update propagation diagrams in 
0(log 2 \) time. 

Edges of G e joining pairs of Steiner points on bisectors are naturally partitioned into 
groups corresponding to ordered triples (b, b 1; f), where b and b x are neighboring bisectors 
with respect to the face f (see Section H] for the definitions). The edges of the initial tetra- 
hedralization T> are assumed to belong to the bisectors incident to them. So, the group of 
edges corresponding to an ordered triple (b, bi, f) consists of all edges from a node on b to a 
node on bi that cross f . Recall that the nodes (Steiner points) on any bisector b were placed 
on a set of segments parallel to the edge of T> incident to b. In our discussion below, we 
refer to these segments, including the edge of V, as Steiner segments. We further partition 

Note that we do not need a priority queue based on elaborated data structures such as Fibonacci heaps. 
Any priority queue with logarithmic time per operation suffices. 
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Figure 11: Two Steiner segments £ and £\ lying on neighboring bisectors b = AABP and 
b x = AACPi respectively, that share a face f = AABC are illustrated. Steiner segments 
on b and bi are parallel to the shared face f . The edges joining nodes on £ and £\ form the 
group of edges corresponding to the triple (£, 

the group of edges associated with the triple (b, bi, f) into subgroups corresponding to pairs 
of Steiner segments (£, £\) on b and b 1; respectively, see Figure [TT] (a). In this way, the 
edges of G e are partitioned into groups corresponding to ordered triples (£,£i,f), where £ 
and l\ are Steiner segments parallel to f on two neighboring bisectors sharing f . The group 
corresponding to (£,£i,f) is denoted by E(£,£ 1 ,f) and consists of all oriented edges from a 
node on £ to a node on £\ that cross f . 

A fixed bisector b has either three or six neighboring bisectors (b itself and two or five 
others, respectively) with respect to each of the two faces forming the dihedral angle bisected 
by b. Hence, the total number of ordered triples (b, bi, f) does not exceed 36n. By Lemma 
13. H the number of Steiner segments on any bisector is 0(-^ log |) and thus the number of 

subgroups of a group corresponding to a triple (b, b 1 , f) is 0{- log 2 -). In total, the number 
of groups E{£, t x , f) is 0(f log 2 ±). 

A node u lying on a Steiner segment £ will have a set of representatives for each group 
E(£,£i,i) corresponding to a triple, where £ is the first element and where its propagation 
set is non-empty. We denote this set by p(u, £i,f). The set of representatives p(u, £\, f) will 
correspond to the structure of the propagation set I(u; £, £\, f), as we will detail in the next 
subsection. 

Consider an iteration of our algorithm. Let v be the node extracted from priority queue 
R, containing all representatives. Let T(v) be the set of triples (£,£i,f) such that v lies on 
£. First, we need to move v from S a to S, since the distance from v to the source s has been 
found. Nodes that are targets of the edges originating from v need to be added to S a . Then, 
we need to compute representatives of v for each group of edges, where edges originating at 
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v participate and where its propagation set is non-empty. Finally, we need to compute new 
representatives for all nodes in the set of predecessors of v, which we denote by R~ x (v ). The 
outline of our algorithm is as follows. 



2. Insert v in S and update S a ; 

3. For each triple (£,£i,f) £ T{v) do 

3.1 Update the data structures related to the Propagation Diagram X(£,£i,f); 

3.2 Find new representatives for nodes whose propagation set has changed in 3.1; 

3.3 Update sets of representatives p(u;£',£,F), for all u £ 

3.4 Update R with respect to 3.2 and 3.3. 

In the remainder of this section, we address the implementation of this algorithm and analyze 
its complexity. First, we observe that the number of iterations is \V £ \. The total number of 
representatives cannot exceed the number of oriented edges in G £ , which is less than \V £ \ 2 
and so, the size of the priority queue R is bounded by | V £ \ 2 (later we show that it is actually 
0(^7=- log 2 |)). Therefore, a single priority queue operation takes 0(log|V^|) time and the 

total time for Step 1 is 0(\V e \ log ]V^|). The total time for Step 2 is 0(|T4| log J). 

In Section I5.2.2[ we describe the structure and maintenance of the data structures re- 
lated to the propagation diagrams X(£,£i,f). Computation and updates of the sets of rep- 
resentatives are described in Section 15.2.31 We conclude our discussion in Section 15.2.41 by 
summarizing the time complexity of the algorithm and by establishing our main result. 

5.2.2 Implementation of Step 3.1 

We consider a fixed triple (£,£±,f), where £ and l\ are Steiner segments on neighboring 
bisectors b and bi sharing f. The propagation diagram X(£,£i,f), was defined as the set 
consisting of the propagation sets of the active nodes on I. Instead of explicitly computing 
the propagation diagram, we construct and maintain a number of data structures that allow 
efficient computation and updates of representatives. 

Consider an iteration of our algorithm. Denote the currently active nodes on £ by 
Ui, . . . ,Uk, and assume that they are listed by their order of entering S. We denote this 
set by S(£) and assume that it is stored and maintained as a doubly linked list ordered 
according to the position of the nodes on I. In Step 3.1, we update the data structures 
related to the propagation diagram X(£,£i, /). According to our definition, the propagation 
set I{u) = I(u;£,£i, f) of a node u £ t consists of all edges (it, «i) in E(£,£i,f) such that 
5(u) + c(it, vi) < S(ui) + c(wj, Vi), for i = 1, . . . , k. Clearly, the set I(u) can be viewed and 
described as a subset of the set of nodes v\ on l\ that satisfy the following three conditions: 

CI. The nodes u and V\ are adjacent in G £ by an edge that crosses f; 

C2. 8{u) + c(it, v i) < S(ui) + c(ui, Vi), for i — 1, . . . , k; 

C3. The node V\ is in S a . 

We construct and maintain separate data structures for the nodes on l\ satisfying each of 
these three conditions: The data structure related to CI is called Adjacency Diagram and is 
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denoted by A(£, £%, /). It consists of sets A(u, £i), for all nodes u on £, where the set A(u, £\) 
consists of the nodes on £\ that satisfy CI. This data structure is static. The data structure 
related to C2 is, in fact, a dynamic additive Voronoi diagram on £\ for the active nodes on 
£ with respect to the weighted distance function c(u, x) defined and studied in Section [2j 
see ([1]). Finally, the nodes on i\ that are in S a are stored in a dynamic doubly-linked list 
and organized in a binary search tree with respect to their position on l\. We denote this 
data structure by S a {£\). The lists S(£) and S a {£\) are readily maintained throughout the 
algorithm in logarithmic time per operation. Next, we describe in detail the construction 
and maintenance of these data structures. 

Adjacency Diagram: The Adjacency Diagram A(£,£\,f) consists of sets A(u,£\), for all 
nodes u on £. We assume that the nodes on £\ are stored in an ordered list V[£\) according 
to their position on that segment. For any fixed node u G £, the adjacency set A(u,£i) will 
be computed and stored as a sublist of the list V(£i). We denote this sublist by A(u,£i). 

We reduce the size of A(u,£i) by replacing each portion of consecutive nodes in them 
by a pair of pointers to the first and to the last node in that portion. (Isolated nodes are 
treated as portions of length one.) Hence, each sublist A(u,£\) is an ordered list of pairs 
of pointers identifying portions of consecutive nodes in the underlying list V(£i). The size 
of the sublists implemented in this way is proportional to the number of the consecutive 
portions they contain. Next, we discuss the structure of the lists A(u,£\) and show that 
their size is bounded by a small constant. 

According to our definitions (Section H]), an edge (u,u\) is present in A(u,£i) if the local 
shortest path if(u,ui] f) does not touch the boundary of /, where the path 7r(u,ui; /) was 
defined in ( I22j) . We refer to intervals on £\ with both of their end-points being Steiner points 
as Steiner intervals. Furthermore, we say that a Steiner interval is covered by the set A(u, £\) 
if all Steiner points, including its end-points, are in A(u,£i). Clearly, each maximal interval 
covered by A(u,£i) corresponds to and defines a portion of consecutive nodes on £\ that are 
adjacent to u. Moreover, by our definition, the list A(u, t\) consists of the pairs of pointers 
to the end-points of the maximal intervals covered by A(u, £i). In the next lemma, we show 
that there are at most seven maximal Steiner intervals covered by A(u,£i). 

Lemma 5.1 The number of the maximal intervals covered by A(u,£i) is at most seven. 
The corresponding ordered list A(u,£\) can be computed in 0(logK(£i)) time, where K[£\) 
denotes the number of Steiner points on £\ . 

Proof: Presented in Appendix 3. □ 

We assume that the nodes that are end-points of the maximal Steiner intervals covered 
by the sets A(u, £±), for all nodes u G £\, are pre-computed in a preprocessing step and stored 
in the lists A(u, £\) as discussed above. Lemma I5TT1 implies that this preprocessing related to 
the group (£,£\,f) takes 0(K(£)\ogK(£i)) time, where K(£) and K(£i) denote the number 
of the nodes on £ and £\, respectively. Next, we discuss the Voronoi diagram data structure 
related to condition C2. 

Dynamic Additive Voronoi Diagram: We assumed that the currently active nodes, 
iti, . . . ,Uk on £, are listed by order of their insertion into S. So, for the distances of these 
nodes to the source, we have S(ui) < • • • < S(uk). We view the distance 8{uj) as an additive 
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weight assigned to the node u i: and consider the additive Voronoi diagram of ui, . . . ,Uk on 
£i with respect to the weighted distance function, introduced and studied in Section [2] and 
defined by flTJ. From the definition (see ([TJ), the weighted distance c(u, x) for a node u on 
i and a point x G i\ is given by 

c(u, x) = c(u, x; f) = min {w|ua| + Wf\aa\\ + w\\aix\}, 

a,a\dF 

where F is the plane containing the face f; w, W\ are the weights of the cells containing I 
and £i, respectively, and Wf is the weight associated to the face /. An important observation 
for our discussion here is that if a; is a node on t\ adjacent to u, then the cost of the edge 
(u, x) is c(u, x). 

We denote the end-points of the segment l\ by A\ and B\ and assume that it is oriented 
so that A\ < B\. For i — 1, . . . , k, the Voronoi cell V(itj) is defined as the set of points on i\ 

V(ui) — {x G (Ax, Bi) : 5(ui) + c(ui,x) < 5(uj) + c(uj,x) for j ^ i}, 

where ties are resolved in favor of the node that has entered S earlier. Clearly, the Voronoi 
diagram V(ui, . . . , Uk) is a partitioning of (Ai, B\) into a set of intervals, where each interval 
belongs to exactly one of the Voronoi cells. Hence, V(u\, . . . , Uk) is completely described by 
a set of points A\ = x < X\ < ■ ■ ■ < x m < x m+ \ = B\ and an assignment between the 
intervals (xj, Xj + i), for j = 0, . . . , m, and the cells of the diagram. 

We assume that V(wi, . . . , Uk) is known and stored. We further assume that a node v on 
I has been extracted by the extract-min operation in Step 1 of our algorithm. In Step 3.1, 
we need to add the new site v and to compute the Voronoi diagram V(u\, . . . ,Uk,v). Next 
we show how this can be achieved in 0(log 2 -) time. First, the following lemma shows that 
the Voronoi cell of v has a simple structure. 

Lemma 5.2 Let ui, . . . ,Uk be the active nodes on i and let v be the last node inserted in 
S. Then the Voronoi cell V(v), in the Voronoi diagram V(u\, . . . ,Uk,v), is either empty or 
consists of a single interval on l\ . 

Proof: By our assumptions 5(ui) < S(v), for i = l,...,k. The Voronoi cell V(v) can 
be represented as an intersection V(v) = Df =1 Vj(u), where the sets Vi(v) are defined by 
Vi(v) = {x G t\ : 5(v) — 5{ui) + c(v,x) < c(ui,x)}. By Theorem 12.11 each of Vi(v) is either 
empty or is an interval on £ 1; and thus the same is true for their intersection. □ 

Using the above lemma, we easily obtain a bound on the size of the Voronoi diagrams. 

Corollary 5.1 The number of the intervals comprising the diagram V(u\, . . . , Uk) does not 
exceed 2k — 1 . 

Next, we present and analyze an efficient procedure which, given the Voronoi dia- 
gram V(ui, . . . , Uk) and a new node v inserted in S, determines the Voronoi diagram 
V(ui, . . . ,Uk,v). This includes computation of the Voronoi cell V(v), update of the set 
of points x±, . . . , x m describing V(u\, . . . , Uk) to another set describing V(iti, . . . , Uk, v ) and 
update of the assignment information between intervals and Voronoi cells. 
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A\ = xq x\ x x'i x + X3 x\ = B\ 




u\ 113 v u-2 



Figure 12: The figure illustrates updates of the diagram V. The Voronoi diagram V(ui, u 2 , u 3 ) 
for the nodes u\, u 2 , and u 3 is characterized by the sequence {xo < X\ < x 2 < 23 < £4} 
and the assignment V(ui) = (xq,xi) U (0:3, X4), V{u 2 ) = (x 2 ,xs), V(u 3 ) = (xi,x 2 ). After 
computation of the Voronoi cell V(t>) = (x",x + ) the Voronoi diagram V(«i, u 2 , 113, v ) is 
characterized by the sequence {xq < X\ < x~ < x + < X3 < X4} and the assignment 
V(«i) = (x ,xi) U (x 3 ,x 4 ), V(u 2 ) = (x + ,x 3 ), V(u 3 ) = V(v) = (x~,x + ). 

According to Lemma l5\2l the Voronoi cell V(v ) is an interval, which we denote by (x~ , x + ). 
Let M be any of the points x\, . . . , x m characterizing the diagram V(ui, . . . , Uk). The fol- 
lowing claim shows that the relative position of M with respect to the interval (x~,x + ) can 
be determined in constant time. 

Claim 5.1 The relative position of M with respect to the interval (x~ ,x + ) can be determined 
in 0(1) time. 

Proof: By the definition of point M, it follows that there are two nodes and Ui 2 such 
that ^(/UjJ + c(ui t , M) = 5{u i2 ) + c(-Uj 2 , M). We denote the latter value by d(M) and note 
that d(M) < 5(ui) + c(ui,M), for i — 1, . . . , k. Then, we compute the value d(v,M) = 
5(v) + c(v, M) and compare it with d(M). 

If d(v,M) < d(M), then we have M e (x~,x + ) and thus x~ < M < x + . In the case 
where d(v,M) > d(M), we compute the Voronoi cell A(v) of v in the three cites diagram 
V(ui t , Ui 2 , v). By Lemma [5.2[ the cell A(t>) is an interval on t\. Since M must be outside 
A(-y) and (x~,x + ) C A(-y), it follows that the relative position between M and (x~,x + ) is 
the same as the relative position between M and A(i>). 

The claimed time bound follows from the described procedure, which besides the constant 
number of simple computations, involves a constant number of evaluations of the function 
c(-, •) and eventually solving of the equations c(-Uj, x) — c(v, x) = 5(v) — S(ui), for i = i 1 ,i 2 .n 

We derive the following binary search procedure, which computes the Voronoi cell V(v). 
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ALGORITHM: Voronoi cell V(v) 

Input: The sequence X = {A\ = xq < X\ < ■ ■ ■ < x m < x m+ \ - 
Output: Points x~ and x + such that V(v) = (x~,x + ). 

A. Compute the point x~ first by the following: 

1. While \X\ > 2 do Steps 1.1 - 1.3 below 

1.1. Find the median M of the sequence X. 

1.2. Determine the relative position between M and x~ . 
1.3 If x~ < M then set X = {x < ■ ■ ■ < M} else set X 

2. If \X\ = 2 compute x" directly. 

B. Compute the point x + in the same way. 

Once the cell V{v) = (x~ ,x + ) has been computed, the update of diagram V(ui, . . . ,Uk) 
to diagram V(«i, . . . ,Uk,v) can be done in a natural way. The sorted sequence of points 
X(ui, . . . ,Uk,v) characterizing the diagram V(ui, . . . ,Uk,v) is obtained from the sequence 
X(ux, . . . ,Uk) = {%o < ■ ■ ■ < x m+ i} by inserting the points x~ and x + at their positions 
and by deleting points (if any) x~ < Xj-+i < • •■ < < x + lying inside the interval 

(x~,x + ), where Xj- and the left and the right neighbors of the points x and x + , 

respectively. We need to delete each of the intervals (xj,Xj + i), for j = j~, . . . , j + — 1, from 
the cell that contains it and to add intervals (xj-,x~) and (x + ,Xj+) to the cells that have 
contained intervals (xj- , and (xj+-i,Xj+), respectively. Indeed, if the cell of some of 

the nodes u±, . . . ,Uk becomes empty then this node is removed from the set of active nodes 
in the group E(£,£i, f). 

To implement all of these updates efficiently, we maintain the sequence X of points 
characterizing the Voronoi diagram of the currently active nodes on £ in an order-statistics 
tree, allowing us to report order statistics as well as insertions and deletions in 0(log|X|) 
time. Based on this data structure, computation of the interval (x~,x + ) takes 0(log 2 \X\) 
time, since it takes 0(log|X|) iterations, and each iteration takes 0(log|X|) time. The 
update of the Voronoi diagram requires two insertions and j + — j~ + 1 deletions in X, where 
insertions take 0(log |X|) time and deletions are done in amortized 0(1) time. 

Let us estimate the time for the maintenance of the Voronoi diagram of the active nodes 
in the group E(£,£±,f). We denoted the total number of the nodes on £ by K{£). Each of 
the nodes on £ becomes active once during the execution. Thus, each node on £ becomes 
subject of the procedure Voronoi cell exactly once. According to Corollary I5.1| the sizes 
of the sequences X characterizing Voronoi diagrams in the group E(£,£i,f) are bounded 
by 2K(£) + 1. Therefore, the total time spent by the procedure Voronoi cell in the group 
E(£,£ u f) is 0{K{£)\og 2 K{£)). In total, there are 0{K{£)) insertions in the sequence X, 
and the total number of deletions, clearly, is at most the number of insertions. Hence, the 
total time spent for insertions and deletions is 0(K(£)\ogK(£)). Thus, the time spent for 
the maintenance of the Voronoi diagram in a fixed group E(£,£i,i) is 0(K(£) log 2 K{£)). 
Next, we discuss the computation and maintenance of a data structure that combines the 
adjacency diagram and the Voronoi diagram. 

Propagation Diagram: As discussed above, the propagation set I(u) of an active node u 
on £ is described completely by the set of nodes on £\ satisfying conditions CI, C2, and C3. 
We denote the set of nodes on £\ satisfying CI and C2 with respect to u by I'(u). Slightly 
abusing our terminology, we refer to this set again as propagation set of u. Similarly, we refer 



Bi}. 



= {M <■■■ < x m+1 }. 
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to the set consisting of the sets I'(u) for all currently active nodes as propagation diagram 
and denote it by X'(u\, . . . , Uk), where, as above, U\, . . . ,Uk are the currently active nodes on 
h- 

The difference between the originally defined propagation set I(u) and the set I'(u) is 
that the elements of I(u) are the edges joining u to the nodes on £\ satisfying CI, C2, and 
C3, whereas the elements of I'(u) are the nodes on £\ that satisfy CI and C2, but not 
necessarily C3. Indeed, the set I'(u) is closely related to I(u) and when combined with the 
list S a (£i) describes it completely. Based on this observation, we compute and maintain the 
propagation diagram X'(u\, . . . , u^j instead of the originally defined diagram. 

We describe the sets I'(ui) by specifying the maximal Steiner intervals they cover. We 
implement these sets as ordered lists of pairs of pointers to the end-points of these intervals 
in the underlying list V(£\). The propagation sets of different active nodes do not inter- 
sect, and hence, the end-points of the maximal Steiner intervals of the propagation sets 
I'(ui), . . . ,I'(uk) form a sequence, X! = {Ai < y\ < z\ < ■ ■ ■ < y mi < z mi < B{\, where 
l\ = (Ai,Bi). The points yj and Zj, for j = 1, . . . ,mi, are Steiner points (nodes) on £\. 
Any of the Steiner intervals (yj, Zj) is a maximal Steiner interval covered by one of the sets 
I'(ui), . . . , I'(uk), whereas the Steiner points inside the intervals (zj, j/j+i) do not belong to 
any of the sets I'(ui). Clearly, the sequence X' plus the assignment of the intervals (yj, Zj) to 
the sets X'(ui) covering them determine the diagram X'(u\, . . . , We implement sequence 
X! as an ordered list of pointers to the underlying list V(£\). In addition, we associate with 
it a binary search tree based on the position of the Steiner points on the segment £\. The 
diagram X'(u\, . . . , itfe) is maintained in Step 3.1 and details are as follows. 

Let, as above, v be the node extracted by the extract-min operation in Step 1 in the 
current iteration of the algorithm. We assume that the diagram X'(u\, . . . , Uk) is known 
- i.e., we know the sequence X' as well as the assignment of the intervals (yj,Zj) to the 
propagation sets I'(uj). Next, we describe the update of X' and the assignment information 
specifying X'(u\, . . . ,Uk,v). By definition, I'(v) consists of the nodes on £\ that lie in the 
Voronoi cell V(v) and belong to the adjacency set A(v,£i). By Lemma [5.21 V(v) is either 
empty or a single interval, which we have denoted by (x~,x + ). We denote by (v~,v + ), the 
largest Steiner interval inside the interval (x~,x + ). The interval (v~,v + ) is easily found 
using binary search in 0(logK(£i)) time, where as above K(£i) denotes the number of 
Steiner points on £\. On the other hand (see Lemma |5TT|) . the adjacency set A(v, £\) consists 
of the nodes lying inside constant number (at most seven) of Steiner intervals, which were 
computed and stored as the list A(v,£i). Hence, the maximal Steiner intervals specifying 
the propagation set I'(v) can be obtained as the intersection of intervals in A(v,£\) with 
(v~,v + ). This is done in constant time by identifying the position of the points v~ and v + 
with respect to the elements of the list A(v,£\). Clearly, the so-computed maximal Steiner 
intervals covered by I'(v) are at most seven. We update the sequence X! by inserting each of 
the maximal Steiner intervals covered by I'(v) in the same way as we inserted the interval 
(x~,x + ) into the sequence X describing the Voronoi diagram. More precisely, let (y, z) be 
any of the maximal Steiner intervals covered by X'(v). We insert the points y and z at their 
positions in the ordered sequence X' , and then we delete the points of X' between y and z. 
If the interval containing y is (yj, Zj), we set new Zj to be the Steiner point preceding y on 
£\. Similarly, if the interval containing z is (yj,Zj), then we set yj to be the Steiner point 
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following z on £\. 

At each iteration of the algorithm, the endpoints of at most seven intervals are inserted 
into the sequence X'. Hence, the size of the sequence X' is bounded by 14K(£) and insertions 
in X' are implemented in 0(logK(£)) time. Deletions are implemented in 0(1) time. The 
total number of insertions is 0(K(£)) and the total number of deletions is at most the number 
of insertions. Therefore, the total time spent for the maintenance of X' and the propagation 
diagram is O (K (£) (log K(£) + K{i x ))). 

Finally, we summarize our discussion on the implementation of Step 3.1. The compu- 
tations and times related to a fixed triple (£,£i,f) are as follows. First, in a preprocessing 
step the lists A(u,£\), for all nodes u on £, are computed in 0(K(£) \ogK(£\)) time (Lemma 
15. ip . Times spent for the maintenance of the lists S(£) and S a (£i) are 0(K(£) \ogK(£) and 
0(K(£\) logK(£i), respectively. The time spent for maintenance of the Voronoi diagram for 
the active nodes on £ requires 0(K(£) log 2 K(£)) time. The time for the maintenance of the 
Propagation Diagram is 0(K(€){\ogK(£) + \ogK(£i))). Therefore, the total time for the 
implementation of Step 3.1 is 



{0(K(t)(\og 2 K(£) + \ogK(£ 1 )) + 0(K(h) log K^))) 




where we have used Lemma [3.11 to estimate that the number of triples (£,£i,i) with a fixed 
first or second element is O(-^log^), and that \ogK(£) and logi^(£i) are O(log^). 

Lemma 5.3 The total time spent by the algorithm implementing Step 3.1 is 0(^=Tog 3 ~). 
5.2.3 Computation and updates of set of representatives 

Next, we concentrate on the computation of representatives in Steps 3.2, 3.3 and 3.4. The 
set of representatives p(v;£, £i,e) of an active node v on £ in a group E(£, £\,f) contains 
one representative for each interval (yj,Zj) in the propagation set I(v). Recall that I(v) 
consists of a set of intervals (yj, Zj) stored in the sequence X', characterizing the propagation 
diagram of the currently active nodes on £. The representative in p(v; £, £\,e), corresponding 
to (yj, Zj) G I{v), is the target of the minimum cost edge from v to a node in S a n (yj, Zj). By 
Lemma [2. 11 the function c(v,x) is convex and thus in any interval it has a single minimum. 
Let x*(v) be the point on £\, where c(v,x) achieves its minimum. To efficiently compute 
the representatives, we compute in a preprocessing step the points x*(v), for all nodes on £. 
From the definition of the function c(v,x) and Snell's law, it follows that x*(v) is the point 
on £\ that is closest to v. So, each of x*(v) can be computed in constant time, which leads 
to 0(K(£)) preprocessing time for the group E(£, £\, f), where K(£) is the number of nodes 
on £. Thus, the total time for preprocessing in all groups is 0(^- log -). 
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We have associated two data structures to the set of nodes in S a that lie on a fixed 
Steiner segment i\. First, we maintain them in a doubly-linked list and second, we maintain 
them in a binary-search tree, with respect to their position on l\. We show that finding a 
representative p{v) G p(v;£,£i,e) takes O(log^) time. There are three situations, where we 
need to compute or update p{v): 

1. New representatives p(v) are computed when v becomes active and its propagation set 
is non-empty. We need to compute one new representative for each maximal Steiner 
interval (y, z) in the propagation set I{v). Recall that there are at most seven such 
intervals and they were computed and stored in the sequence X 1 . 

To compute p(v) in the interval (y, z), we determine the leftmost and rightmost nodes 
from S a inside the interval {y,z). This is done by finding the position of the points 
y and z in the sequence of nodes currently in S a . Let the leftmost and the rightmost 
nodes from S a in (y, z) be y a and z a , respectively. Then, we determine the position of 
the point x*(v) with respect to y a and z a . 

If it is to the left of y a , then p(v) = y a . If it is to the right of z a , then p(v) = z a . If 
x*(v) is inside (y a ,z a ), we determine the two nodes in S a immediately to the left and 
to the right of x*(v), and p(v) is one of these two nodes. Using the binary-search tree 
on S a , the nodes y a and z a and eventually the nodes neighboring x*(v ) are determined 
in 0(log i) time. 

2. When some representative p(v) is removed from S a , a new representative for v is one 
of the neighbors of p(v ) in the doubly-linked list S a that lie in the same interval (y, z) 
as p{v). This is done in 0(1) time. 

3. When some interval of the propagation set I(v) shrinks and the current representative 
p(v) is no longer inside this interval, then p(v) is updated as follows. Let, as above, y a 
and z a be the leftmost and the rightmost nodes from S a , respectively, in the updated 
interval. Then, if p{v) lies to the left of y a , we set p(v ) = y a . If p(v ) is to the right 
of z a , we set p(v) = z a . As above, determination of the nodes y a and z a is done in 
O(log^) time. 

To complete our analysis, we need to estimate the total number of representatives which 
are computed by our algorithm. Each pair (representative, predecessor) relates to the edge 
joining them. Since such a pair can be computed at most once by the algorithm, the total 
number of representatives related to nodes that are vertices of T> is bounded by the total 
number of edges incident to these nodes, which is 0(|V^|). It remains to estimate the number 
of representatives which are related to nodes that are Steiner points. Consider an iteration 
for a node v that is a Steiner point. There are O(-^log^) triples in 7~(f), and at most 
nine new representatives are computed in Step 3.2. For each predecessor in i? _1 (w) that is 
a Steiner point, a single representative is computed. The number of predecessors |i? _1 (v)| 
is 0(^| log i). Hence, in a single iteration, O(^log^) representatives related to Steiner 
points are computed. Since the number of iterations is OQV^I) and the computation of a 
single representatives takes 0(log ~) time, we obtain that the total time for the execution of 
Steps 3.2 and 3.3 is 0(H log 2 ±). 
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Finally, the number of priority queue operations executed in Step 3.4 is bounded by the 
number of computed representatives. Thus, the total time for Step 3.4 is 0{jQ log j log ~). 

5.2.4 Complexity of the algorithm and the main result 

Here, we summarize our discussion from the previous three subsections and state our main 
result. Step 1 of our algorithm takes 0(|K|log^) time. Step 2 requires 0(|14|) time. By 
Lemma 15.31 Step 3.1 takes in total 0(^=Tog 3 ^) time. The total time for implementation 
of Steps 3.2 and 3.3 is log 2 ±) and the total time for Step 3.4 is log f log ±). By 

Lemma [3.11 we have that \V e \ = O(^jTog^). We have thus established the following: 

Theorem 5.1 The SSSP problem in the approximation graph G £ can be solved in 
CK^logf log 3 i) time. 

Consider the polyhedral domain T>. Starting from a vertex vq of T>, our algorithm solves 
the SSSP problem in the corresponding graph G e and constructs a shortest paths tree rooted 
at vq. According to Theorem I4.1[ the computed distances from vq to all other vertices of 
T> (and to all Steiner points) are within a factor of 1 + e of the cost of the corresponding 
shortest paths. Using the definition of the edges of G £ , an approximate shortest path can 
be output by simply replacing the edges in the discrete path with the corresponding local 
shortest paths used to define their costs. This can be done in time proportional to the 
number of segments in this path, because computation of the local shortest paths takes 0(1) 
time. The approximate shortest paths tree rooted at Vq and containing all Steiner points 
and vertices of V can be output in 0(|V^|) time. Thus, the algorithm we described solves 
the WSP3D problem and the following theorem states the result. 

Theorem 5.2 Let V be a weighted polyhedral domain consisting of n tetrahedra and e G 
(0,1). The weighted shortest path problem in three dimensions (WSP3D), requiring the 
computation of approximate shortest paths from a source vertex to all other vertices ofT>, 
can be solved in Of-Jg- log - log 3 -) time. 



6 Conclusions 

This paper generalizes the weighted region problem, originally studied in 1991 by Mitchell 
and Papadimitriou [21] for the planar setting, to 3-d weighted domains. We present the 
first polynomial time approximation scheme for the WSP3D problem. The complexity of 
our algorithm is independent of the weights, but depends upon the geometric features of the 
given tetrahedra as stated in Lemma [3.11 

There are some fairly standard techniques which can be employed here to remove the 
dependence on geometry (cf., [1]), provided that there is an estimate known on the maximum 
number of segments (i.e., the combinatorial complexity) in weighted shortest paths in three 
dimensions. It can be shown that the combinatorial complexity of weighted shortest paths 
in planar case is Q(n 2 ) |2T|. We conjecture that the same bound holds in three dimensions, 
but the proof techniques in [21] do not seem to apply here, since they use planarity. If the 
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combinatorial complexity of these paths in three dimensions is a polynomial in n, then we 
can remove the dependence on the geometry by increasing the run time by a polynomial 
factor in n. We do not recommend this approach, since the increase in the running time will 
be significant. Already, in the planar case (and in terrains), in an experimental study |18j . 
it was shown that a constant number of Steiner points suffice to produce high-quality paths. 
We believe that the same holds here and this merits further investigation. 

This paper also investigated additive Voronoi diagrams in heterogeneous media. We 
studied a fairly simple scenario and already the analysis of that was very technical and 
cumbersome. It is desirable to find simpler and more elegant ways to understand the combi- 
natorics of these diagrams. Nevertheless, we believe that the discretization scheme and the 
algorithms presented here can be used successfully for efficient computation of approximate 
Voronoi diagrams in heterogeneous media. 

Our algorithm does not require any complex data structures or primitives and as such 
should be implementable and even practical. Its structure allows Steiner points to be gen- 
erated "on the fly" as the shortest path wavefront propagates though the weighted domain. 
This feature allows the design of more compact and adaptive implementation schemes that 
can be of high practical value. 

One of the classical problems that motivated this study is the unweighted version of this 
problem, namely the ESP3D problem. There, we need to find a shortest path between a 
source and a target point, lying completely in the free space, avoiding three-dimensional 
polyhedral obstacles. We can use our techniques to solve this problem, though this will 
require triangulating (i.e., tetrahedralization) the free space. As outlined above, the com- 
plexity of our algorithm depends upon the geometry of these tetrahedra; so it is natural to 
ask whether the free space can be partitioned into nice tetrahedra? Unfortunately, there is 
no simple answer to this question which has been an important topic of study in computa- 
tional and combinatorial geometry for several decades. Nevertheless, our algorithm provides 
a much simpler and so far the fastest method for solving the ESP3D problem, provided the 
free space is already partitioned into non-degenerate tetrahedra. 

Combining the techniques of answering weighted shortest path queries on polyhedral 
surfaces [2] and the existence of nice separators for well-shaped meshes [20], we believe that 
our construction presented in this paper can be used for answering (approximate) weighted 
shortest path queries in 3-d. 

References 

[1] Pankaj K. Agarwal, R. Sharathkumar, and Hai Yu. Approximate Euclidean shortest 
paths amid convex obstacles. In Claire Mathieu, editor, SODA, pages 283-292. SIAM, 
2009. 

[2] Lyudmil Aleksandrov, Hristo Djidjev, Hua Guo, Anil Maheshwari, Doron Nussbaum, 
and Jorg-Riidiger Sack. Algorithms for approximate shortest path queries on weighted 
polyhedral surfaces. Discrete & Computational Geometry, 44(4):762-801, 2010. 



43 



[3] Lyudmil Aleksandrov, Anil Maheshwari, and Jorg-Riidiger Sack. Approximation algo- 
rithms for geometric shortest path problems. In STOC, pages 286-295, 2000. 

[4] Lyudmil Aleksandrov, Anil Maheshwari, and Jorg-Rudiger Sack. Determining approxi- 
mate shortest paths on weighted polyhedral surfaces. J. ACM, 52(l):25-53, 2005. 

[5] Tetsuo Asano, David Kirkpatrick, and Chee Yap. Pseudo approximation algorithms 
with applications to optimal motion planning. Discrete & Computational Geometry, 
31(1):139-171, 2004. 

[6] Franz Aurenhammer and Rolf Klein. Handbook of Computational Geometry, chapter 
Voronoi Diagrams. North Holland, 2000. 

[7] Chandrajit L. Bajaj. The algebraic degree of geometric optimization problems. Discrete 
& Computational Geometry, 3:177-191, 1988. 

[8] John F. Canny and John H. Reif. New lower bound techniques for robot motion planning 
problems. In FOCS, pages 49-60. IEEE, 1987. 

[9] Joonsoo Choi, Jurgen Sellen, and Chee-Keng Yap. Approximate Euclidean shortest 
paths in 3-space. Int. J. Comput. Geometry Appl, 7(4):271-295, 1997. 

[10] Joonsoo Choi, Jurgen Sellen, and Chee-Keng Yap. Precision-sensitive Euclidean shortest 
path in 3-space. SI AM J. Comput, 29(5):1577-1595, 2000. 

[11] Kenneth L. Clarkson. Approximation algorithms for shortest path motion planning 
(extended abstract). In STOC, pages 56-65. ACM, 1987. 

[12] B. Cox and B. Treeby. Artifact trapping during time reversal photoacoustic imaging for 
acoustically heterogeneous media. IEEE Trasaction on Medical Imaging, 29(2):387-396, 
2010. 

[13] Sariel Har-Peled. Constructing approximate shortest path maps in three dimensions. 
SIAM J. Comput, 28(4):1182-1197, 1999. 

[14] Philip L. Inderwiesen and Tien- When Lo. Fundamentals of Seismic Tomography, vol- 
ume 6 of Geophysical Monograph Series. 1994. 

[15] Rolf Klein. Concrete and Abstract Voronoi Diagrams, volume 400 of Lecture Notes in 
Computer Science. Springer, 1989. 

[16] Rolf Klein, Elmar Langetepe, and Zahra Nilforoushan. Abstract Voronoi diagrams 
revisited. Computational Geometry, 42(9) :885 - 902, 2009. 

[17] J. Krozel, S. Penny, J. Prete, and J.S.B. Mitchell. Comparison of algorithms for synthe- 
sizing weather avoidance routes in transition airspace. Collection of Technical Papers - 
AIAA Guidance, Navigation, and Control Conference, 1:446-461, 2004. 

[18] Mark Lanthier, Anil Maheshwari, and Jorg-Riidiger Sack. Approximating shortest paths 
on weighted polyhedral surfaces. Algorithmica, 30(4):527-562, 2001. 



44 



[19] Ngoc-Minh Le. Abstract Voronoi diagram in 3-space. Journal of Computer and System 
Sciences, 68(1):41 - 79, 2004. 

[20] Gary L. Miller, Shang-Hua Teng, William Thurston, and Stephen A. Vavasis. Automatic 
mesh partitioning. In Alan George, John Gilbert, and Joseph Liu, editors, Graphs 
Theory and Sparse Matrix Computation, The IMA Volumes in Mathematics and its 
Application, pages 57-84. Springer- Verlag, 1993. Vol 56. 

[21] Joseph S. B. Mitchell and Christos H. Papadimitriou. The weighted region problem: 
Finding shortest paths through a weighted planar subdivision. J. ACM, 38(l):18-73, 
1991. 

[22] Joseph S. B. Mitchell and Micha Sharir. New results on shortest paths in three dimen- 
sions. In SoCG, pages 124-133. ACM, 2004. 

[23] Christos H. Papadimitriou. An algorithm for shortest-path motion in three dimensions. 
Inf. Process. Lett, 20(5):259-263, 1985. 

[24] J.H. Reif and Z. Sun. Movement planning in the presence of flows. Algorithmica (New 
York), 39(2): 127-153, 2004. 

[25] Zheng Sun and John Reif. Bushwhack: An approximation algorithm for minimal paths 
through pseudo-Euclidean spaces. In Peter Eades and Tadao Takaoka, editors, Algo- 
rithms and Computation, volume 2223 of Lecture Notes in Computer Science, pages 
160-171. Springer Berlin / Heidelberg, 2001. 10.1007/3-540-45678-3_15. 

[26] Zheng Sun and John H. Reif. On finding approximate optimal paths in weighted regions. 
Journal of Algorithms, 58(1):1 - 32, 2006. 

[27] T. Varslot and G. Taralsen. Computer simulation of forward wave propagation in 
soft tissue. IEEE Trasactions on Ultrasonics, Ferroelectrics, and Frequency Control, 
52(9):1473-1482, 2005. 



15 



A Appendix 



A.l Proof of Proposition 12.21 

Proposition ET^t The second mixed derivative of the function x = x(y,a) is negative, i.e. 

Xya ^ • 

Proof: First, we consider the case where w~ = w + . In this case, the function x(y, a) can 
be represented and differentiated explicitly. Recall, that the path 7f(t>,x) in this case either 
consists of a single segment or is a three segment path as shown in Figure [3] (b). So, in 
the case where the path consists of a single segment, we have x(y, a) = y cot a. In the case 
where the path consists of three segments, x(y, a) = y cos a/y/n 2 — cos 2 a, where n = w/w~ . 
The mixed derivatives x ya of these two functions are —1/ sin 2 a and — k 2 sin a/ (k 2 — cos 2 a) 2 , 
respectively, and both are readily negative. 

Next, we consider the case where w~ 7^ w + . We introduce some additional notation as 
necessary for our presentation below (Figure [2]). We denote the coordinates of the bending 
point a of the path tt(v, x) by a = (x~ , y + ). Furthermore, we set x + = x—x~ and y~ = y—y + . 
Clearly, x~ , x + , y~, and y + can be viewed as functions of the independent variables y and 
a. We have x = x~(y~(y, a), a) + x + (y + (y, a), a) and thereby 

x y = X y-Vy + x t+yy- ( 3? ) 
We differentiate f[3"T|) with respect to a and obtain x ya = A + A\ + A 2 , where 

A = X y-Vm + x t+yt^ A i = x y- y -y*yy + x y+ y +vtv^ and A 2 = x~. a y~ + x+ +a y+. 

We complete the proof by showing that the terms A, A±, and A2, are negative. 

We begin with the term A = x~_y~ a + Xy + y ya . From the identity y = y~ + y + , it follows 
that y~ a + y ya = and hence, A = — x^ + )y ya . From our notation, Snell's law, and the 
relation cos a = cos 9 simp, we derive the following: 

, z + cos a , z + v sin 2 tp — cos 2 a 

x = . — ; y — = — , 

\J k 2 — sin 2 if yK 2 - sin 2 <p 

and (38) 

z~ cos a z~ \J sin 2 09 — cos 2 a 

x = y = — = — . 

V 1 — sin 2 - sin 2 </? 

First, we compute x~_ and We denote sin 2 ip by a, differentiate x + and y + with respect 
to a and obtain x+ + as the ratio We differentiate expressions (138]) with respect to a 

and obtain 

r „ < / x ^ , 1 m 4- z + cosa , 2 + (k 2 — cos 2 a) . , 

{2 G "H : c(«, x) + C < c(t/, = — - T , y+ = 1 (39) 

— <t)2 ivo" — cos^ a{K z — a) •■ 
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and hence 



, , . , cos ol\J sin 2 ip — cos 2 a , i<v . 

x+ = a£/y+ = . 40 

y k 1 — cos 2 a 



16 



Similarly, 



z cos a _ z 1 - cos a 

x„ = 5-, y n = — — o-, and 

2(1 -<j)i 2vV-cos 2 a(l-a)i 

cos a a/ sin 2 ip — cos 2 a 

V = = 1 2 ■ 41 

y 1 — cos"' a 



So, for the difference x y _ — x y+ we have 



x — x, = cos ay sin 2 ip — cos 2 a[ 



1 1 



y y v i _ cos 2 a k 2 — cos 2 a ' 



_ , ,2 n cosay/sinV - cos 2 a 
(1 — cos a){K' L — cos z a) 

The latter shows that 

sign(ar_ - = sign(/€ 2 - 1). (43) 

To prove the negativity of A = (x~_ — x y+ )y~ a , we show that sign(y~ Q ) = sign(l — k 2 ). We 
have 

Vy = ValVo = Va/iVa + Vt) = , , + / _ and thuS Vya = ~ 7 ' T7 ^2 ' 

Hence, it is sufficient to show that 

sign((y+/y;) a ) = sign(K 2 - 1). (44) 

So, we continue by establishing the sign of the derivative (y£ /y~) a . We use fl39|) and fj4Tj) 
and compute the ratio y^/y^ as follows 

+ / - ^ + (^ 2 -cos 2 a)(l -<x) 3/2 , + / 

= —q 2 w 2 ^72" = A )BC, where 

2; (1 — cos z a){K z — <j) 6 ' z 

n K 2 -cos 2 a fl-a\ i/2 ( 1 - sin 2 ip \ 3/2 

5 = — and C = = — rV 1 • (45) 

1 — cos z a \k z — a J \hi z — sin p> J 

Then, we compute the derivatives B a and C a using the expressions (H5|) 

g sin 2cc(l — cos 2 a) — sin 2«(k 2 — cos 2 a) (1 — k 2 ) sin 2a 
(1 — cos 2 a) 2 (1 — cos 2 a) 2 

3 / 1 — sin 2 <p \ ^ 2 (— sin 2(p) (k 2 — sin 2 9?) — (— sin 2ip)(l — sin 2 tp) 
a 2 \k 2 — sin 2 <p J (k 2 — sin 2 ip) 2 ^ a 

3 / 1 - sin 2 </? \ 1/2 sin 2y?(l - k 2 ) 



2\k 2 — sin 2 ip J (k 2 — cos 2 a) 2 
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and obtain 



[z-/z + ){y+/y-) a = B a C + BC a (46) 

3 1 

1 — sin 2 tp \ * (1 — k 2 ) sin 2a 3 ( 1 — sin 2 tp ^ 5 (k 2 — cos 2 a) sin 2</>(l — n 2 )tp a 



k 2 — sin 2 tp J (1 — cos 2 a) 2 2 \ k 2 — sin 2 (/? / (1 — cos 2 a)(/t 2 — sin 2 tp) 2 

= (1 - ;: 2 ) Vl-sin 2 ^ ^ 

(k 2 — sin 2 tp) 3 / 2 {\ — cos 2 a) ' 

sin 2a(l — sin 2 </>) 3 sin2y9(K 2 — cos 2 a) 

where £> = 1 — r^— tp a . 

1 — cos 2 a 2(k 2 — sin tp) 

Omitting the positive multiplicative terms in (l4"6"j) . we derive that sign(y£ /y~) a = sign((l — 
k 2 )D) and continue with the evaluation of sign(D). We compute the derivative tp a using the 
identity y = y~ + y + , which implies = y~ + y+. We differentiate expressions from (1381) 
with respect to a and obtain 

(1 — sin 2 <p) sin 2a + sin 2<p(l — cos 2 a)ip a 

y a = I 2 )- 



y + a = {z + m 



(sin 2 tp — cos 2 a) 1 / 2 (l — sin 2 tp) 3 / 2 
[k 2 — sin 2 tp) sin 2a + sin 2tp(n 2 — cos 2 a)<^ c 



(47) 



(sin tp — cos 2 a) 1 / 2 (/« 2 — sin tp) 3 / 2 
From these two, we obtain 

/ sin 2a (1 — sin 2 (p)(k 2 — sin 2 (p) 

tp a = ^ , ; A ^, where (4* 

J sin 2tp 



I = Z -( K 2 -sm 2 tp) 1/2 + 2 + (l -sinV) 1/2 and (49) 
J = z~(l- cos 2 a)(K 2 - sinV) 3/2 + ^ + (k 2 - cos 2 a)(l - sinV) 3/2 . 

Next, we substitute tp a from (I48p in the expression D given in (1461) and obtain 

_ . /, • 2 > 2 J — 31 (k 2 — cos 2 a)(l — cos 2 a) 

D = smacosafl — sin tp) — -. 

J (1 — cos 2 a) 

The term sinacosa(l — sin 2 tp) and the denominator in this expression are positive and by 
expanding the numerator we have 

sign(D) = sign[2J — 3I(k 2 — cos 2 a)(l — cos 2 a)] 

= sign [2z~(l - cos 2 a)(n 2 - sin 2 tp) 3 / 2 + 2z + (k 2 - cos 2 a)(l - sin 2 tp) 3 / 2 

— 3z~(l — cos 2 a)(n 2 — cos 2 a)(n 2 — sin 2 tp) 1 / 2 — 3z + (l — cos 2 a)(n 2 — cos 2 a)(l — sin 2 y?) 1 ^ 2 ] 

= sign [z~(l - cos 2 a)(/t 2 - sin 2 tp) 1/2 D~ + z + (k 2 - cos 2 a)(l - sin 2 <p) 1/2 D + ] , 

where D~ = 3 cos a — 2 sin tp — k and D = 3 cos a — 2 sin tp — 1. (50) 
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Now, we observe that the terms multiplied by D + and by D~ are positive and show that 
D~ and D + are negative. We use cos a = cos # sin and cosa K = cos#sin</? re , where sin<y9 = 
k sin ip K and obtain 

I 

D = 3 cos a — 2 sin 99 — 1 = 2 cos a — 2 sin ip — sin a 

O/-1O O O OO^i o 

= 2 cos b'sin — 2 sin (p — sin a = — 2 sin <£>sin 6* — sin a<0 

O 000/0 0\ 

and D~ = 3 cos a — 2 sin tp — n — k (3 cos a K — 2 sin y> K — 1) = 

k 2 (— 2 sin 2 (/? K sin 2 — sin 2 a K ) < 0. 

From (130]) and -D + , Z) _ < O, we get sign(D) = — 1. From f|46|) it follows that sign ) Q = 
— sign(l — k 2 ) and thus siga(y~ ) = sign(l — « 2 ). The latter implies that A < 0. 

Next, we consider the term A±. From the identity y~ + y+ = y Q = 0, we get Ai = 
VaiXy-y-Vy ~ x y+ y +Vy)- To evaluate the sign of y~, we substitute ip a from (1481) in the 
expression (147|) and by omitting positive multiplicative term, we obtain 

sign(y~) = sign[( J — (1 — cos 2 a)(«; 2 — sin 2 <p)I) cos a] (51) 
= sign{z + (l — sin 2 ^^[(k 2 — cos 2 a)(l — sin 2 tp) — (1 — cos 2 a)(ft 2 — sin 2 cos a} 
= sign[(l — k 2 ) cos a]. 

Next, we evaluate the sign of the difference x~_ y _y y — x y+y+ y y . We compute x y+y+ as follows 

+ , + . . + cos a . z + (k 2 — cos 2 a) 

X y + y + = ^ X v+' a /y^ = 7rr~o o \ / = t = / 



2(k 2 — cos 2 ^)"^ — cos 2 a 2\J o — cos 2 «(k 2 — a) 2 



where we have differentiated (l4"0]) with respect to cr = sin 2 <£> and used ( 1391) . We compute 

_ in the same way and obtain 

y y J 

+ cos a(K 2 — sin 2 <^)i _ cos q(1 — sin 2 tp)? 

% y + y + = z+{^- C o^af and Vr = z -(l-cos 2 a) 2 " (52) 

Furthermore, we have y+ = y+/y a = y+/ (y+ + y~) and y" = y+/ (y+ + y~). So, we compute 
y+ and y~ using (1551) as follows 

, z + (k 2 — cos 2 a)(l — sin 2 </?)§ _ — cos 2 a)(/t 2 — sin 2 <p)% . . 
< = j and y y = , (53) 

where J was defined in (j4"9~j) . Using (1521) and (1531) . we determine 

sign(x y _ y _y y -x+ +y+ y+) = sign[(l/(l -cos 2 a) - 1/(k 2 -cos 2 a)) cos a] = sign[(ft 2 - 1) cos a}. 

The latter and ( IBTj) imply A\ < 0. 

Finally, we show that A 2 = x~^ a y y + x y+a y y is negative too. We first compute the 
derivative x y+a by differentiating the expression ( 1401) with respect to a. We have 



+ P K sin a + cos a sin 9? cos9?(k 2 — cos 2 a)<p a 

X y + a 1-1 9 \i/ 9 9 \9 ^ ' 

(sin <£> — cos z a) 2 («r — cos z a) z 
49 



where P K = 2 cos 2 ai^n 2 — sin 2 ip) — sin 2 ip(n 2 — cos 2 a). We substitute ip a using the expression 
( l4"gj) and multiply by y+ using the expression (|5T>|) . After simplification, we obtain 

, , , . . 9 ,o/n JP K — cos 2 a(n 2 — cos 2 a) (I — sin 2 </?)(k 2 — sin 2 ip) I 

x + y + = z + sin a 1 - sin 2 w 3/2 — i J -— r ^ 55 

J^sm — cos^ a) 2 (k z — cos z a) 

Analogously, we obtain the following 

. 9 9 , o/ 9 JPi — cos 2 ail — cos 2 — sin 2 ip)(n 2 — sin 2 ip)I . , 

x _ y„= z sma(K -smV ' — i — — , (56) 

j [sm ip — cos a) 2 (1 — cos a) 

where Pi = 2cos 2 a(l — sin 2 ip) — sin 2 y9(l — cos 2 a). We sum (1551) and (1561) . simplify and 
omit the positive multiplicative terms obtaining 

sign(x"_ a ?/- + = (57) 

sign[^(K 2 - sin 2 v?) 3/2 (K 2 - cos 2 a)Qi + z + (l - sinV) 3/2 (l - cos 2 a)Q K ], 

where 

Qi = J P\ — cos 2 a(l — cos 2 a) (1 — sin 2 ip) (k 2 — sin 2 ip)I and 
Q K = JP K — cos 2 a(i% 2 — cos 2 a)(l — sin 2 (p)(K 2 — sin 2 ip)I . 

We denote the expression in the square brackets by R. Finally, evaluate R and show that 
it is negative. First, we evaluate and simplify Q K and Q\ . We substitute the expressions / 
and J from (1491) in Q K and group the terms containing z~ and z + . Then, we substitute the 
expression for P R from (l54"j) and by simplification we get 

Q K = z~(k 2 — sin 2 ip) 3 / 2 [(l — cos 2 a)Pk — cos 2 a(n 2 — cos 2 a)(l — sin 2 ip)] + 

z + (l — sin 2 p) 3 / 2 [(k 2 — cos 2 a)Pfc — cos 2 a(n 2 — cos 2 o)(k 2 — sin 2 <^)] 
= z~(k 2 — sin 2 <y9) 3//2 (cos 2 a — sin 2 p)(n 2 + cos 2 a — 2k 2 cos 2 a) + 
z + (l — sin 2 </?) 3//2 (cos 2 a — sin 2 ip)k 2 (k 2 — cos 2 a). 

In the same way, we obtain the following representation for Qi 

Qi = z~(k 2 — sin 2 <y9) 3 / 2 (cos 2 a — sin 2 tp)(l — cos 2 a) + 

z + (l — sin 2 <p) 3 ^ 2 (cos 2 a — sin 2 <p)(n 2 + n 2 cos 2 a — 2 cos 2 a). 

Substitution of Q K and Q\ in (|5"T|) produces an expression for R of the form 

R = (z^) 2 R 1 + (z + ) 2 R 2 + z~z + R 3 , where 
Pi = (k 2 — sin 2 (p) 3 (n 2 — cos 2 a)(l — cos 2 a) (cos 2 a — sin 2 ip) 
P 2 = (1 — sin 2 </j) 3 (l — cos 2 a)n 2 {n 2 — cos 2 a)(cos 2 a — sin 2 99) 
P 3 = (1 - sin 2 <^) 3/2 (k 2 - sin 2 <^) 3/2 (cos 2 a - sin 2 ip) x 

[{k —cos a)(/t +k cos a — 2cos a) + (1 — cos a){n + cos a — 2k cos a)J 
= (1 - sin 2 ipj 3/2 (K 2 - sin 2 (^) 3/2 (cos 2 a - sin 2 y?) x 

[(k 2 - cos 2 a) 2 + k 2 (1 - cos 2 a) 2 + (1 - k 2 ) 2 cos 2 a] 

From these expressions, it is clear that terms Pi, P2, and P3 are negative, since cos 2 a — 
sin 2 ip < and all other terms are positive. Thus, P and consequently A 2 are negative. The 
proposition is proved. □ 
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A. 2 Upper bound on the constant Cabp(^) 
Next, we show that 

11IABI , A\AB\ 2 h 

C A Bp(t) < / x ■ 9/ /o\ lQ g2 



r(e)sin 2 ( 7 /2) 62 r(e)r(A)r(B) ' 

Recall that A = (1 + sin( 7 /2)) -1 and < e < 1. We use the following two inequalities, 
which are easily derived from the properties of logarithms: 
For any X > 1 

3.44 lnX , X 2, 

log A - 1 X< , ln-<ln-log 2 X 58 

Vesin(7/2) £ £ 

By our definition of the constant CABpif) an d ( IT7|) it follows that 

2e|AB| |AB| e|AB| 

CAMt) -r(e)X(l- A) log f ^ gA T^FTO + r(e)(l - A) 2 log f (59) 

2£ 2 , /i 4£ 2 , \AB\ 

+ \ 2 lo gA-i— TT + i 2 lo £A- 



logf er(e) logf &A ~ eAv/r(A)r(B)' 

We estimate the terms on the right-hand side of this inequality using inequalities ( 158]) above. 
For the first one, we have 

2e\AB\ \AB\ < (2^2 + l) 2 y/e\AB\ 2\AB\ 

r(e)A(l - A) logf 0Sa_1 £ \ y/r(A)r(B) ~ V2r{e) sin( 7 /2) log f ° gA_1 e^/r(A)r(£) 

3.44(2^+ 1) 2 |AB| , 2IABI \AB\ , 2|AB| , s 

< F v Q — — tt In — t= < 25 ; ' 9 ' r logo— ^ ■ 60 

~ v / 2r(e)sin 2 ( 7 /2)logf e^IgB) r(e) sin 2 ( 7 /2) S2 ^HA)HB) 

For the second one, we have 

e\AB\ < (2y^+l) 2 |^| [ABj 
r(e)(l - A) 2 logf _ r(e)sin 2 ( 7 /2)logf r(e) sin 2 ( 7 /2) log f ' 1 J 

The sum of the third and fourth terms is estimated by 



2£ 2 /. h \AB\ \ Ue 1 - 5 f, I h . 2\AB 



log A -i— r T + 21og A -i— = < — — — 2 InW— T + ln 



logf y er(e) eA^r^r^) / sm( 7 /2)logf\^ \j er(e) e y /r(A)r(B) i 

14£ L5 ln2 1 2\AB\Vh £ L5 , A\AB\ 2 h , , 

< ; — logo — , < 5 ; — -—log 9 —- — -— r — ; — -. (62) 

- sin( 7 /2) 62 yJr{e)r{A)r{B) sin( 7 /2) S2 r(e)r{A)r{B) V ; 

We substitute Q5QJ), (EI]), and fl52) in (EH]), use 2r(e) < |AB|, r(e) < h and obtain 

\AB\ , 4|AB| 2 /i 
CUbp(*) < 23 '. ' . . log ■ 



r(e)sin 2 ( 7 /2) ta r(e)r(A)r(B) ' 
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A. 3 Proof of Lemma 15.11 



Lemma 15.11 The number of the maximal intervals covered by A(u,£i) is at most seven. The 
corresponding list A(u,£\) is computed in 0(logK(£i)), where K(£i) denotes the number of 
Steiner points on l\ . 

Proof: First consider the case when the segments i and £\ lie in different tetrahedra. We 
denote the weights of the tetrahedra containing £ and £\ by w~ and w + , respectively. Let F 
be the plane defined by the face / and let F~ and F + be the two half-spaces defined by F, 
where we assume that £ is in F~ and £\ is in F + . Furthermore, we assign weights w~ and 
w + to F~ and F + , respectively, and we consider shortest weighted path 7t(u,y) between u, 
that is on £, and an arbitrary point y on l\. 

As discussed in Section[2l in this case, the path tt(u, y) has the form 7f (u, y) = {u, a(y),y}, 
where the point a(y) lies in F and is uniquely defined by Snell's law (Figure [3] (a)). By our 
definition, there is an edge joining u to a Steiner point U\ G £\ if and only if the point a{u\) 
lies in the interior of the triangle /. The interior can be represented as intersection of three 
half-planes defined by the lines containing the sides of /. So, we first obtain an upper bound 
on the number of maximal intervals covered by A(u,£\) in the case where / is a half-plane 
defined by an arbitrary line L in F. 

There is one-to-one correspondence between end-points of the maximal intervals and the 
points y on £\ for which a(y) lies on L. Hence, the number of maximal intervals covered 
by A(u,£i) can be estimated by counting the number of points y for which a(y) lies on L. 
When the point y traverses the segment £i, the bending point a(y) defines a curve in the 
plane F, which we denote by a(y). We consider Cartesian coordinate system MiV in F, such 
that segment £\ projects onto the segment (^i, [12) of the //-axis, and u projects onto y-axis, 
say, at the point u' = (0,u ). We denote the /i-coordinate of the projection of y by fi(y) 
or simply by fi when no ambiguity arises. Then, as we discussed in Section [2J the curve 
a.(y) has a representation a(/x(y)) = (rfi(y), (1 — t)u ), where r is the unique solution of the 
equation (j2J). 

Let L(fi, v) be the linear function, such that v) — represents the line L in 0^ v . 
Then, a point a(y) belongs to L if L(r/x(?/), (1 — t)vq) — 0. Thus, we obtain the following 
system of algebraic equations for r and ji{y) 



where z + is the Euclidean distance from £\ to F and z~ is the Euclidean distance from 
u to F. Excluding r from this system leads to a degree four algebraic equation for fi(y). 
Therefore, there may be no more than four intersections between a(y) and L and, hence, the 
number of the maximal intervals covered by A(u,£i) in the case where / is a half-plane is at 
most three. 

In the case where / is a triangle, we denote by f\, f 2 , and ^3 the three half-planes defining 
/ and by I\, I2, and ^3 the corresponding sets of maximal intervals. Then, any maximal 
interval A defined by / is obtained as an intersection Aifl A 2 fl A3, where Aj £ Jj, % = 1, 2, 3. 
Each of the sets Jj contains at most three intervals and it is easily seen that the number of 
intervals that are intersections of the type Ai n A 2 H A 3 does not exceed seven. 




(63) 
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-Un 1*0 /' / M 



Figure 13: The figure illustrates the curves a(y) and ai(?/) in the case (a) when v~ + v + > u 
and in the case (b) when v~ + v + < Uq. A line L and its intersections with the curves a(y) 
and B.i(y) are shown. The number of maximal intervals covered by E(v,£i) in the case when 
/ is the lower half-plane defined by L, is 2 for both instances (a) and (b). 



The list A(u,£i) is computed by finding the position of the solutions of the systems (1631 
with respect to the elements of V(£i). This can be done by performing binary search and 
hence the computation of the list A(u,£i) in this case takes 0(logK(£i)) time. 

Next, we consider the case where u and £\ lie in the same tetrahedron, say the one that is 
in F~ . If the weight w~ is smaller than w + , then the path n(u, y) has the form {u, a(y),y}. 
The point a(y) is the point in F that lies on the segment (u,y'), where y' is the point 
symmetric to y with respect to F. So, the curve a(y), in this case, is a segment. Hence, 
each of the sets /j, for i — 1, 2, 3, in this case, is either empty or consists of a single interval. 
Consequently, there can be at most one interval obtained as intersection of intervals in these 
sets. Thus, in this case, there can be no more than one maximal interval covered by A(u, £±). 

Finally, we consider the case where w~ is greater than w + . In this case, the shortest 
path 7t(u,y) has the form {u,a(y),ai(y),y}, where the segment (a(y),ai(y)) lies in F. We 
discussed the structure of this path in Section [2] and illustrated it in (Figure [3] (b)). The 
curves a(y) and &i (y) have explicit representations as we detail below. We set v~ = z~ tan if* 
and z/ + = z + tany9*, where the critical angle ip* is defined by siny?* = w + /w~ and z~ , z + 
are the distances from v and £\ to F, respectively. 

In this notation the curves a(y) and a\(y) have the following representations in 0^ v 

a(y) = { " = ^ fOT ^ l</i0 (64) 

u = v - yjjy ) 2 - /J 2 for /Jo > \fi\ < v , 

3Ll (y) = { V= ^f^ fOT l^l<^0 ^ (65) 

(U ~ U)y / (U+) 2 -V 2 = flU for \fj,\ > /I Q , 

where fio = z - +z+ yj [y~ + v + ) 2 — if v~ + v~ > u and fi = if V + v~ < u . 

In the case where v~ + v~ > v (Figure [TBI (b)), the curve a(y) is a half-circle centered 
at the point (0, u ) with radius v~ . The curve ai(y) is symmetric with respect to the z/-axis. 
The part of ai (y) to the right of O v is monotonically decreasing. It is convex and approaches 
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the /^-axis at infinity. In the case where v~ + v~ < u (Figure [13] (a)), the curves a(y) and 
SLi(y) have a common part - a horizontal segment that projects at (— /io,/io) on the /i-axis. 
For > /i , the curves are the same as in the case u~ + V > u . 

Again, we consider, first, the case when / is a half-plane defined by a line L in F. There 
is an edge between u and a point y on l\ if and only if the segment (a(y), ai(y)) lies entirely 
inside /, i.e. in one of the half-planes defined by L. By a simple case analysis, we determine 
that in the case where / is a half-plane, the number of maximal intervals covered by A(u, i\) 
is at most 2. 

Hence, in the general case when / is a triangle, each of the sets I\, I2, and ^3, as defined 
above, contains at most 2 intervals. The number of intervals that can be intersections of 
the type flf =1 Aj with Aj £ Jj is at most 4. The latter proves that the number of maximal 
intervals covered by A(u,£i) in the case where I and l\ are in the same tetrahedron, whose 
weight w~ is bigger than the weight w + of the neighboring tetrahedron, is at most 4. 

Computation of the list A(u,£i) in this case is done again by binary search and takes 
0(K(£i)) time. The lemma is proved. □ 
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